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Abstract 

In the preceding paper, linear response methods have been applied to obtain formally exact 
expressions for the parameters of Navier-Stokes order hydrodynamics. The analysis there is general, 
applying to both normal and granular fluids with a wide range of collision rules. Those results are 
specialized here to the case of smooth, inelastic hard spheres with constant coefficient of normal 
restitution, for further elaboration. Explicit expressions for the cooling rate, pressure, and the 
transport coefficients are given and compared with the corresponding expressions for a system 
of elastic hard spheres. The scope of the results for further analytical explorations and possible 
numerical evaluation is discussed. 
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I. INTRODUCTION 



Granular materials frequently exhibit flows similar to those of normal fluids, and for prac- 
tical purposes these flows are often described by phenomenological hydrodynamic equations 
[1, 2}. In the simplest cases, such equations have the form of Navier-Stokes hydrodynamics 
with an energy sink representing collisional "cooling". The parameters of these equations 
(cooling rate, pressure, and transport coefficients) are unknown in general. For normal fluids, 
methods of non-equilibrium statistical mechanics, e.g. linear response, have been applied to 
derive formally exact expressions for these parameters in a form that is suitable for introduc- 
ing approximations. In the preceding paper [3] , this analysis has been extended to granular 
fluids with similar results. The objective of the present work is to make those results more 
explicit by specializing to the idealized model granular fluid: a system of smooth, inelastic 
hard spheres or disks. As has been shown in ref. [3], this idealization is a limiting form 
for more realistic collisions, where the interaction energy scale, e.g. energy of deformation, 
is small compared to characteristic kinetic energies. This can be controlled by sufficient 
activation of the fluid, so the predictions of this model are more than academic [4]. Also, as 
in the case of simple atomic fluids, this idealized model of hard particles provides the most 
tractable setting for further theoretical analysis of the results obtained in this work. 

A significant advantage of this limit is a scaling of the reference homogeneous cooling 
state (HCS) to which linear response is applied. This scaling property results in several 
simplifications to the formalism described in [3]. For example, a stationary representation 
for the statistical mechanics of this granular fluid follows from the use of dimensionless phase 
space variables, where the particle velocities are scaled with the cooling thermal velocity in 
the HCS ensemble [5, 6]. Furthermore, the hydrodynamic response functions for small spatial 
perturbations of the HCS are given in terms of stationary time correlation functions in this 
representation. These time correlation functions are composed from two biorthogonal sets 
of observables, referred to here as the direct and conjugate functions. The direct functions 
are the microscopic observables whose ensemble averages are the hydrodynamic fields. The 
conjugate functions are generated by functional derivatives of a local HCS with respect 
to the hydrodynamic fields. For hard spheres or disks, it is possible to transform these 
functional derivatives into ordinary phase space derivatives. Generally, the simplifications 
made possible for hard spheres allow a somewhat more transparent interpretation of the 
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dynamics and the structure involved in the formal expressions for transport coefficients. 
They also provide potentially better access to evaluation of these expressions using molecular 
dynamics (MD) simulations. 

The purpose of this companion paper is twofold. The first is to specialize the results 
obtained in [3] to the particular system of inelastic hard particles, to give explicit expressions 
for the various hydrodynamic parameters, and to take advantage of this simplified setting 
to interpret their content [7]. The second is to record the various tools of nonequilibrium 
statistical mechanics of a system of inelastic hard particles whose utility transcends the 
context of their application in this work. Regarding the first purpose, details of the analysis 
are suppressed in the main text, with the reader referred to the preceding paper or to 
the Appendices here. Instead, emphasis is placed on the final results, and similarities or 
differences between Helfand and Green-Kubo representations for normal and granular fluid 
transport coefficients. 

The second purpose is accomplished in most of the eight Appendices, which summarize 
the statistical mechanics and dynamics for inelastic hard spheres. Although much of the 
formalism appears in part elsewhere, there are many new results as well. In addition to the 
generators for dynamics of phase space observables and the Liouville equation, the generator 
for time reversed dynamics is identified. Also, the microscopic conservation laws for both 
the forward and time reversed dynamics are given, and the two sets of different fluxes are 
identified. Although straightforward to obtain, these results do not appear in the literature 
even for elastic hard spheres. Similarly, although the Helfand representation for transport 
coefficients has been used in MD simulations of elastic hard spheres, the corresponding 
Green-Kubo expressions have not been discussed explicitly until recently [8, 9]. All such 
elastic hard sphere coefficients are given in detail in Appendix F using the more general 
granular fluid formalism developed here. Dimensionless forms of the Liouville equation, 
including one suitable for MD simulation [10, 11], the microscopic conservation laws, and a 
new set of balance equations for the conjugate functions noted above are derived in these 
Appendices. All of these results are critical for interpreting the exact expressions given 
here for the transport coefficients, but are also the building blocks for more general future 
investigations of granular fluid properties using the non-equilibrium statistical mechanics of 
hard spheres. 

The layout of this presentation is as follows. In the next section, the primary results of 
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the previous paper are specialized for the case of inelastic hard particles whose collisional 
loss in energy is characterized by a constant coefficient of normal restitution a. Only a broad 
outline of the results is given in the main text, with the details of the dynamics and statistical 
mechanics of hard spheres recorded in the Appendices. Then, explicit expressions for the 
various hydrodynamic parameters are given. The cooling rate and pressure are obtained as 
averages over the HCS, which also can be expressed in terms of integrals over the associated 
reduced two particle distribution function. The transport coefficients are given both Helfand 
[12] and Green-Kubo [13] representations, in terms of the long time limit of appropriate 
stationary time correlation functions for the HCS. Next, the relationship of these results to 
those obtained from kinetic theories is sketched briefly, as well as the possibility of numerical 
evaluation of the results using molecular dynamics (MD) simulation. Finally the key points 
in this work are summarized and some concluding remarks are made. 

II. LINEAR RESPONSE 

The microscopic model considered here is a system of smooth mono-disperse hard spheres 
(d = 3) or disks (d = 2) of mass m and diameter a, whose collisional loss in energy is 
characterized by a constant coefficient of normal restitution a. The details of the model, 
the statistical mechanics of this system, the generators of dynamics, and time correlation 
functions are given in Appendix A. In this section, the simplifications to the linear response 
theory developed in the previous paper, due to this choice of collision model, are identified 
and the results there translated to the case at hand. 

A. Homogeneous Cooling State 

The reference state with respect to which the linear response of this system is studied, is 
the homogeneous hydrodynamic state whose entire time dependence arises due to the cooling 
temperature. This is the so-called homogeneous cooling state (HCS). First this reference 
state is characterized in the context of both the phenomenological macroscopic hydrody- 
namics, and the more fundamental statistical mechanics. At the level of hydrodynamics, 
this state is described by a solution with constant number density n h and flow velocity Uh, 
and a homogeneous but time dependent temperature T h (t), that obeys the hydrodynamic 



4 



equation 

{dt + Co[n h ,T h (t)]}T h (t) = 0, (1) 

where (o[rih,Th(t)} is a cooling rate that must be specified from the microscopic theory. 
By means of a Galilean transformation, it is always possible to consider Uh = 0. For the 
case of inelastic hard spheres or disks considered here, there is no microscopic energy scale 
associated with the collision model. Therefore, the temperature dependence of the cooling 
rate can be determined by dimensional arguments as Co [nh,T h (t)} oc T^ 2 (t). Equation (1) 
then can be integrated to obtain the time dependence of the temperature, 
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being a dimensionless cooling rate and a "thermal velocity", respectively. Moreover, / is 
an appropriate length scale; for example, the mean free path. Equation (2) is the familiar 
Haff's cooling law [2] that is well established for inelastic hard sphere fluids as one of the 
signatures of the HCS. 

At the level of statistical mechanics, the reference ensemble corresponding to this 
macrostate, the HCS ensemble, is given by a homogeneous "normal" solution to the Liouville 
equation, i.e., one whose entire time dependence occurs through the cooling temperature. 
Furthermore, the absence of any additional microscopic energy scale for hard spheres, im- 
plies that this temperature dependence can occur only through the scaling of the particle 
velocities v r with the thermal velocity v [T h (£)], 

-Nd * I I 9rs v r ~~ Uf t 



Ph [T;T h (t)}) = [lv (t)])- Nd pl 
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where T = {q r , v r ; r — 1, . . . , iV} is a point in the phase space of the system and q rs = 
q r — q s the relative position of particle r with respect to particle s. This special form of 
the iV-particle distribution function allows the temperature dependence of many average 
properties such as the cooling rate, pressure, and transport coefficients, to be determined 
without explicit calculation. Another interesting consequence of the scaling nature of the 
HCS ensemble is its restriction to a constant total momentum surface in phase space, 
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Here S(x) is the Dirac delta function and P = J2 r mv r the total momentum. The proof of 
this result is given in Appendix B. 

It is useful to represent the dynamics more generally for other states in terms of 
the same dimensionless scaled variables as occur in the HCS, i.e. in terms of T* = 
{q*, v*; r — 1, . . . , N}, with q* = q r /l and v* = v r /v (t), where v (t) is defined from the 
temperature of a reference HCS state having the same intial total energy as the system under 
consideration. In this way, it is seen that the HCS is a stationary solution to the dimen- 
sionless Liouville equation, and the HCS time correlation functions also become stationary. 
Thus, this representation for the statistical mechanics, allows the problem of linear response 
about a non-equilibrium time dependent reference state to be mapped onto one where the 
reference is stationary, in closer correspondence to the equilibrium case for normal fluids. 
More details of this transformation are given in Appendix B. 



B. Hydrodynamic Response 

In this subsection, the response of the fluid to small spatial perturbations about the 
reference state described above is considered. The macroscopic variables of interest are the 
dimensionless deviations, 5y* a) of the hydrodynamic fields, y a , from their HCS values, y a ,h- 
They are defined by 

{Sy*J = {Sn*,ST\SU*}, (6) 

5n*- Sn - n ~ nh 5T*~ ST - T ~ Th ^ sir - SU _ u ~ u h m 
" n h n h ' T h (t) T h (t) ' v (t) v Q (t) ' [) 

where n(r,t), T(r,t), and U(r,t) denote respectively the number density, temperature, 
and flow velocity fields. The dynamics of these hydrodynamic fields is given by a response 
equation of the form 

Sy:(k*,s) = J2c:,(k*,s)6y;(k*,0). (8) 

P 

A Fourier representation in the reduce space variable has been introduced to take advantage 
of the homogeneity of the reference state. Also, the dimensionless time scale s defined 
through 

ds = (9) 

is used. The matrix of response functions C*, gives the effect of weak spatial perturbations 
of the hydrodynamic fields on the dynamics of this fluid. Such response functions can 
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be identified in two ways [3]: either from the linearized phenomenological hydrodynamic 
equations or from non-equilibrium statistical mechanics. Equating the two in their common 
domain of validity (large space and time scales), allows identification of exact expressions 
for the parameters of the phenomenological hydrodynamic equations. 

The solution to the dimensionless, linearized phenomenological hydrodynamic equations 
gives 

C*W(k*,s) = e- sK * hvd V°*\ (10) 

If the components of the flow velocity are chosen to be the longitudinal component relative to 
k*, 5C/|* = k* ■ 5U*/k* and d — 1 orthogonal transverse components, 5U* Li) the dimensionless 
transport matrix JC* hyd turns is block diagonal, 

. . / lC hyd \ 

K *hyd II / U N 

^ K? yd ) 

Here ]C* hyd is the "longitudinal" part, corresponding to j<5n*, ST*, 5U*\ j, and it is given by 
/ -ik* \ 

^hyd (r) = Q d^_ + (if _ C «) fc*2 f + (2A1 - (*T) k *2 + ^ k * 

ipj dlnp h ,* ip* h ,* Co , [ 2(d-l) * , *1 7*2 

(12) 

while /C^^ is associated to the "transverse" components, SU±i, which are decoupled from 
the former, 

K,* 2 hyd (k*) = + V *k*^j I, (13) 

where / is the unit matrix of dimension d — 1. The dimensionless parameters and transport 
coefficients in the above equations are defined by 

Ph = -^r> Co = — , A 



/ 



n h Th d ' " ln h Vo 
A* = 7^ — > ^ = ;— > K 



/T^o' mn h lvQ mn h lvo 

c u = c u , r = ^, c t = ThC_ (14) 

/do /D 

Above, Dfe is the pressure of the granular fluid, Co is the cooling rate, A is the thermal 
conductivity, \i is the new transport coefficient associated with granular fluids that measures 
the contribution of density gradients to the heat flux, i] is the shear viscosity, and k is the 



bulk viscosity. Finally, ( u , ( n , and ( T are transport coefficients arising from the local cooling 
rate. The phenomenological definition of all these quantities has been given in ref. [3]. Upon 
introducing the associated dimensionless quantities in Eq. (14), advantage has been taken 
of the already mentioned fact that the temperature dependence of all the parameters in the 
hydrodynamic equations can be determined by dimensional arguments, as a consequence of 
being dealing with hard particles. 

The dimensionless transport matrix JC* hyd is a constant independent of time s and, conse- 
quently, the response function C* is related to the transport matrix by the simple exponential 
form given in Eq. (10). This is one of the primary simplifications that occur for this hard 
sphere or disk system. Further, the eigenvalues of the JC* hyd matrix give the hydrodynamic 
modes of the system, and its eigenfunctions identify the particular excitations of the hy- 
drodynamic fields that give rise to these modes. These can be used to get physical insight 
into the hydrodynamic response of this fluid. Some of these considerations are addressed in 
Appendix C. 

Alternatively, an exact response equation of the form of Eq. (8), can be identified by 
carrying out a linear response analysis at the level of statistical mechanics of the system. In 
this way, the response function C* is found to be a matrix of time correlation functions of 
the form 

C* p (Jfe*; s) = V*- 1 J dT*a* a (r*; k*) e~ sT % (P\ -k*) . (15) 

Here, V* is the volume of the system in reduced units and the "direct" functions a* are 
dimensionless Fourier transforms of linear combinations of the microscopic densities whose 
ensemble average gives the hydrodynamic fields, 

a* a (k*) = r d J dre ikr a* a (r) = J dr* e ik * ^ a* a (r) , (16) 

K} = {^(V-^),e*}, (i7) 

where 

(V ,r l(P}s mWM}. (18) 

I n h n h T h (t) n h mv {t) J 
The phase functions Af, £, and Q are the microscopic number density, energy density and 
momentum density, respectively. Their mathematical expressions are given in ref. [3]. The 
"conjugate" functions ipp in Eq. (15) are generated by functional derivatives of the local 
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HCS ensemble pih [T\ {y a }] as also described in [3], 



ij* (T*;k*) = M p J dre lh - r 



Spin [T\{y a }] 



(19) 

{y a }={n h ,T h (t),0} 



ha (r) 

In the above expression, the Mg's are factors that render the functions ipp dimensionless, 
namely 

M l = [lv (T h )] dN n h , M 2 = [lv (T h )] dN T h) M„ = M ±l = [lv (T h )] dN v (T h ) . (20) 

The local HCS inherits some of the scaling nature of the true HCS ensemble, and is identified 
as [14] 

Pih (r| { Va }) = |n [lvo(q r )r d \ Pi ({^p, Vr ~^ r) } |n) , (21) 

where v (q r ) = v [T(q r )]. The functional dependence on the density field is obtained by 
considering the HCS in an external inhomogeneous field and inverting the relation giving 
the density profile [14]. Clearly, for uniform hydrodynamic fields this becomes the HCS of 
Eq. (4). 

The generator of the dynamics C* in Eq. (15), is a composition of the dimensionless 
Liouville operator L* that generates the trajectories in phase space plus a velocity scaling 
operator, 

r=l r 

The Liouville operator L associated with this system of hard particles is identified explicitly 
in Appendix A. The second operator on the right hand side of Eq. (22) rescales the velocities 
with an s dependence associated with the thermal velocity, and represents the effects of 
cooling in the reference HCS. This simplified form of the modified generator of dynamics in 
the time correlation functions is another special feature of the hard sphere collision model, 
or equivalently, of the absence of any other microscopic energy scale. 

Using the form of the local HCS given in Eq. (21), the functional derivatives with respect 
to the temperature and flow velocity in Eq.(19), can be simplified to the forms 

J s =l L 0± ^ > {y a }={n,T,0} 



J s=1 I ^||,-LU</sj J {j, a }={ n ,T,0} 

AT o 

= -l^^^—pir). (25) 

If this linear response analysis were done for a normal fluid using the canonical Gibbs 
ensemble p c (X), the functions ip* a would reduce to quantities proportional to the familiar 
microscopic densities a Q multiplied by the Gibbs equilibrium ensemble p c . Hence, the dif- 
ferent forms of the conjugate densities appearing in Eq. (15) for the case of granular fluids, 
are a manifestation of the fact that the linear response is now being measured about a 
non-equilibrium macrostate. 

Two different representations for the hydrodynamic response matrix C* have been iden- 
tified in Eqs. (10) and (15), respectively. The latter representation is exact and valid for 
all times and for all length scales. On the other hand, the former representation is obtained 
from the phenomenological hydrodynamic equations, which are the relevant description of 
the fluid on length scales long compared to the mean free path and time scales long com- 
pared to the mean free time. In this limit, the exact response functions (15) must go over 
to that given in Eq.(10), thereby providing an exact identification of all the hydrodynamic 
parameters in terms of time correlation functions. To facilitate this identification, the ex- 
act analog of the hydrodynamic transport matrix K.* hyd defined in Eqs. (11)- (13) above is 
identified as 

/C* (Jfe*, s) = - \d s C* (jfe*; s)l C*- 1 (Jfe*; s) , (26) 
and the cross over to hydrodynamics stated above is 

)C* h y d (k*)= lim lC*(k*,s). (27) 

s»0,fc*«l 

This is the prescription carried out in ref. [3] to obtain Helfand and Green-Kubo represen- 
tations for the transport coefficients of a granular fluid. To lowest order, k* = 0, the exact 
expression for the cooling rate of the fluid is identified. At Euler order, linear order in k*, 
the pressure and the Euler transport coefficient ( u are identified. Finally, to order k* 2 the 
six Navier-Stokes transport coefficients are all obtained. 
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III. HYDRODYNAMIC PARAMETERS AND TRANSPORT COEFFICIENTS 



A. Overview 

In this Section, the exact explicit expressions for each of the parameters occurring in 
the phenomenological hydrodynamic equations at Navier-Stokes order are given, and their 
technical and physical content is noted. The dimensionless units are used throughout, and 
for simplicity the length scale is now chosen by the condition rihl d = 1. 

The required "equations of state" are expressions for the cooling rate and the pressure 
as functions of the density and temperature. These are given as averages of specific phase 
functions over the HCS distribution function. The phase functions are sums of functions 
involving the phase coordinates of pairs of particles, and so the averages can be expressed 
as integrals over the associated reduced two-particle distribution function defined by 



This function depends on the positions of particles only through q{ 2 = q* — q 2 , as a con- 
sequence of the translational invariance of the HCS (homogeneity). The occurrence of 
for the properties considered here, depends on q* 2 only on the sphere for the pair at contact, 
i.e. for q\ 2 = a* = a /I. The velocity dependence can be represented in terms of the relative 
velocity g* 2 = v\ — v 2 , and the center of mass velocity G* 2 = (uj + v 2 ) /2. If the center of 
mass velocity can be integrated out, the result is a probability distribution for the magni- 
tude of the relative velocity and the cosine of the impact angle, i.e., x — \q{ 2 - g^l/ludu a ^ 
contact, 



As shown below, the cooling rate, the pressure, and some parts of the transport coefficients 
can be expressed as low degree moments of F£ (g{ 2l x). 

The dimensionless transport coefficients w* will be given in two equivalent represen- 
tations. The (intermediate) Helfand expressions are long time limits of time correlation 
functions, while the Green-Kubo expressions are given in terms of time integrals of related 
time correlation functions. These two representations are related by the simple identities 




(29) 




(30) 
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with 

n = n H ( s = o), n G (s) = ^h(s). (31) 

The first equality in Eq. (30) identifies the Helfand expression, while the second one is 
the Green-Kubo representation. The symbol lim denotes the ordered limits of V* — > oo 
followed by s — > oo. Note that what is called here the Helfand representation for the sake of 
simplicity, was termed the intermediate Helfand representation in ref. [3]. The correlation 
functions f2#(s) appearing in the Helfand representation have the general form 

n H (s) = v*- 1 J dr* F* s > f (r*) e' s ( r - x ) s m* (r*) . (32) 

The phase functions F* s '- f are defined in terms of either a volume integrated source S or 
a flux / present in the microscopic conservation laws for the direct functions a* defined in 
Eq. (17). The phase functions M*(T*) are first space moments of the functional derivatives 
of the local HCS distribution function. They are obtained from Eq. (19) to first order in k, 
then having the general structure 

'Spm [r| {y a }} 



M*(T*) =M drk-r 



5y(r) 



(33) 

{y a }={n h ,T h (t),0} 



Finally, the dynamics in Eq. (32) is generated by the modified operator C* of the Liouville 
equation, minus one of its eigenvalues \ a determined from 

(r-A Q )c(r*;0)=0. (34) 

As discussed in refs. [3] and [14], the functions ^*(r*;0) are the exact eigenfunctions 
of the modified Liouville operator, representing homogeneous perturbations of the HCS 
hydrodynamics. The eigenvalues \ a = 0,Q/2, —Q/2, the last one being rf-fold degenerate, 
are therefore the same as those of the transport matrix JC* hyd (0). Since the transport 
coefficients characterize the response to spatial variations, it is reasonable that the dynamics 
in their representation as time correlation functions does not include such homogeneous 
excitations. 

The correlation functions f2c( s ) present in the Green-Kubo representation have the cor- 
responding forms 

n G ( s ) = v*- 1 J dr* f* sj (r*) e - s ( r ~ x )r* (r*) , (35) 
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with the new "fluxes" T* (T*) given by 

T* (r*) = - (T - Xj M* (r*) . (36) 

The terminology "flux" for these functions is justified in Appendix E, where it is shown that 
they are the volume integrals of fluxes occurring in the microscopic balance equations for the 
conjugate functions In detail, it is found that the F* s, f (r*) are projected orthogonal 
to the set of eigenfunctions j^* (P*; 0)|, so that the zeros of the operator C* — A do not 
occur in the dynamics. Consequently, the limit on the time integral in Eq. (30) is expected 
to exist. This is the analogue of the "subtracted fluxes" present for the same reason in the 
Green-Kubo expressions for normal fluids [13]. 



B. The cooling rate and the pressure 

The cooling rate was identified in ref. [3] as proportional to the average (negative) rate 
of change of the total energy in the HCS. Its dimensionless value is given by 

C»* = ^7^w.(r> :( r.) (37) 

where 

W*(T*) = - [ dr*w*(T*;r*), (38) 



d _ 

w* being the source term in the dimensionless balance equation for the energy, Eq. (D27). 
Then, using Eqs. (D30) and (D16), it is found that 

1 2 N N 

w (n = J2 £ s( q ; s - a*)e (-<? rs ■ g* rs ) \<? rs ■ g ; s \ 3 . m 

r=l s^r 

Here, cfi. s = q* s /q rs and Q(x) is the Heaviside step function. Therefore, Co is proportional to 
the third moment of the normal component of the relative velocities of all pairs at contact, 
averaged over one of the collision hemispheres, in the HCS, weighted by the fractional loss in 
the kinetic energy of the two particles that is proportional to (1 — a 2 ). Since the cooling rate 
is determined by the average of a sum of phase functions involving only phase coordinates 
of pairs of particles, it can be expressed in the reduced form 

Co = J dq\j dv{J dq* 2 j dv* 2 S(q* 12 - a*)Q(q* 12 ■ g* 12 ) 

x\Q*i2-9l2\ 3 r h {2) (Ql2,<v* 2 ). (40) 
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Transforming to relative and center of mass coordinates, allows further simplification to 

(1 - a 2 )ir d l 2 a* d - 1 f 
Co = 1 J dgl 2 ®{x)g*x*F*{g* l2 , x), (41) 

where x = cos is now the projection of g\ 2 along an arbitrary direction. This rather simple 
result is still exact, and retains all of the two particle (position and velocity) correlations 
representative of the HCS. The two forms (37) and (41) are suitable for complementary 
molecular dynamics simulations of the cooling rate. 

Similarly, the pressure is identified as the HCS average of the trace of the microscopic 
momentum flux 

p* = tL I rfr>*(r*)trH*(r*). (42) 



V*d_ 

Here the phase tensor H* is the volume integral of the dimensionless flux defined in Eq. 
(D28), 



H*(r*) = J dr*h*(T*;r*) 

N (-> , \ * N N 

= E « + 4 E E 5 (& - ^ • 9rs) 

r=l r=l sj^r 

x(« s -^ s ) 2 « s . (43) 
The pressure can be partitioned into two terms, 

P*=Pk+P*c ( 44 ) 

where 

PK = \j dY* p* h {Y*)vf = I (45) 

and 

P*c = (1 2 + ^ / dq{ f dvl J dq* J dv* 5 (q* 12 - a*) (-<? 12 ■ g* u ) 

= ± {d/2) } d J dgl 2 glje(x)x 2 F*(gl 2 ,x). (46) 

The first term on the right hand side of Eq. (44) is the kinetic part of the pressure. It 
arises purely from the transport of momentum associated with the free streaming of the 
particles and implies that px = UhT h (t) . The second term, determined by the two particle 



14 



distribution function at contact, is the collisional transfer part of the pressure. Like the 
cooling rate, it is given by a moment of the simple distribution Fj* (g{ 2 , x). Again, Eqs. (42) 
and (44)- (46) provide different, complementary expressions for the evaluation of the pressure 
by means of molecular dynamics simulations. 

The pressure and the cooling rate, as functions of the hydrodynamic fields n and T, are the 
required equations of state for the hydrodynamics equations which are determined entirely by 
the HCS. This is analogous to the case of normal fluids where the dependence of the pressure 
on n and T is given by the equilibrium Gibbs state. Note that this definition is independent 
of the presence or absence of other normal stresses, which appear as additional independent 
dissipative contributions to the hydrodynamic equations. Similarly, the cooling rate has 
additional independent dissipative contributions characterized by transport coefficients. 



C. Euler Transport Coefficient Q u 

As already indicated, a new transport coefficient, ( u , occurs at Euler order in the tem- 
perature equation of the granular fluid. It is a measure of the contribution to heat transport 
from divergences in the flow velocity, and does not occur in the "perfect fluid" Euler equa- 
tions for a normal fluid that have no dissipation. In Appendix F, the intermediate Helfand 
representation for the dimensionless form of this transport coefficient is identified as 

C^ = hm^( S ), (47) 

where the time correlation function is 

tlx (s) = V*- 1 J dT* W s (r*) e _S ( £ + ^)mIu(T*). (48) 

Here W (F*) is the volume integrated source term W (T*) given in Eq. (39), projected 
orthogonal to the set |^* (r*; 0) j, as discussed above. More precisely, it is 

w s (r*) = (i - p^)w(t*) = w (r*) - ^ (e* - ^n^J - q + n. (49) 

The projection gives the additional terms proportional to the total number of particles iV 



and the total energy E* = ^2 r v* . The moment Ai^u in Eq. (48) is given by 

N a 

^Mn = -E^-^(n, (so) 



dv* 

r=l ' 



15 



that is a space first moment of the velocity derivative of the HCS distribution. To pro- 
vide some interpretation for A4^ Uy note that for a normal fluid, when the canonical Gibbs 
ensemble p c (T) is considered as the reference state, 

N 

Mlu = 2j2q*r-<P*c, (51) 

r=l 

which is the space moment of the momentum density of the fluid. This is modified for 
granular fluids to account for the non-equilibrium nature of the HCS. 

The generator of the dynamics in Eq. (48) consists of the scaled Liouville operator C 
together with the factor Co/2- The latter is one of the eigenvalues (A = — Co /2) solution of Eq. 
(34), and occurs here because M*^j has a non- vanishing component along the corresponding 
eigenfunction in this set of zeros for C* - A. The projected form of W s (T*) assures that 
this component does not contribute to (* u . Further discussion of this effect for the other 
transport coefficients is given below. 

The equivalent Green-Kubo representation is 

(* u = n< u + lim f ds' fif (s') , (52) 
Jo 

with Vtf = fif(O) and 

fig 7 ( s ) = V*- 1 J dT* W s {r*)e~ S ( C + ^) T* (U {r*), (53) 

where 

T*„ = - (T + ^pj M* cu (54) 

is a "conjugate flux", associated with the new conserved quantities of the dynamics of 
this system [3]. The usual integrands of the time integral in Green-Kubo expressions are 
flux-flux correlation functions. Here, one of the fluxes has been replaced by the source in 
the energy equation due to collisional loss of energy. This is characteristic of the Green- 
Kubo expressions for transport coefficients arising from the cooling rate. For all other 
transport coefficients, the familiar flux-flux correlation form occurs, although involving both 
a direct and a conjugate flux. This difference in the two sets of fluxes is due to both the 
dissipation in the collision rule and the singular nature of hard spheres. The occurrence 
of the instantaneous part to the Green-Kubo relation, , has similar origins and is not 
present for normal fluids with non-singular forces. 
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The instantaneous part of the Green-Kubo representation can be simplified further to 

ftf = (V*d)~ 1 J dT*W s (T*)Mlu(T*) = -^(l-a)p* c , (55) 

where p* c is the collisional part of the pressure given in Eq. (46). Hence, (* u is the contribu- 
tion of the source to what would physically constitute the effects of hydrostatic pressure at 
the Euler order, where it enters in the combination (2p*/d) +(* u . If a small volume element 
of the fluid is considered, then the amount of pressure that the fluid element can exert on its 
boundaries is decreased by the energy lost locally due to collisions. Part of the effect of this 
transport coefficient is to decrease the effective pressure in the system. At the level of linear 
hydrodynamics, the two coefficients are indistinguishable in their physical consequences. 



D. Shear Viscosity 

The (intermediate) Helfand and Green-Kubo expressions for the dimensionless shear vis- 
cosity rf are also elaborated in Appendix F to the the form 

rj* = lim n v H (s) = Q v + lim / ds' n v G (s') , (56) 

Jo 

where the correlation functions are defined by 

«h(«) = - d 2 V + * d 1 _ 2 J2J2 J dT * "une'^^M^n, ( 5 ?) 

i=i j=i 

^c(s) = -^^EE / ^H^.(r)e- S ( r+% )TV.(r), (58) 

i=l j=l 

and fig = f^(0). In the above expressions, H*-(r*) is the volume integrated momentum flux 
given by Eq. (43), M-*^ is the traceless tensor 

M;dn = -\ t + '-4 " ~M ■ a|) »* (r *>- (59) 

and the associated Green-Kubo conjugate flux T* ^ is 

T; y = -(r + |^; y .(n. ( 6 o) 

The generator of the dynamics is now i^C* + ^fj and, therefore, has the k = mode of Eq. 
(34) with A = -Co/2 subtracted out. 

17 



In order to better understand the differences from the corresponding expression for a 
normal fluid, note that for the latter with the canonical Gibbs reference state p c (T), it is 

N / o \ 

= J2 U,k> + - -/^r ■ <) pun, m 

r=l ^ ' 

and 

T; >0 -(r) = -TM^in. (62) 

Using the property (A27) and taking into account that L* p* = 0, it is obtained that, for a 
normal fluid, 

TV-(r) = -2 (H*r - I^tr H-) p* c (n, (63) 

where H*~ is the volume integrated momentum flux for the time reversed microscopic con- 
servation laws identified in Appendix D, namely 

N * N N 

H*-(r*) = E^ + ySE^-Oe^-^) 

r=l r=l s^r 

X{<£s ■ grs) 2 qr S qrs- (64) 

Therefore, the normal fluid version of Eq.(56) is 

ot/*— 1 d d „ s „ 



P :(r). 

(65) 

The normal fluid shear viscosity is given by the long time integral of the equilibrium time 
correlation functions of the forward momentum flux H* with the time reversed momentum 
flux H*~, together with an instantaneous part. This instantaneous term remains nonzero in 
the elastic limit, as an artifact of the singular nature of the hard particle interaction. Note 
that in the elastic limit, the hydrodynamic modes defined through Eq. (34) all have the 
same eigenvalue A = 0, so the generator for the dynamics in both the elastic and inelastic 
case can be thought of as having the k* = mode subtracted out. 

The instantaneous contribution Qq for the granular fluid can be expressed in terms of the 
reduced pair distribution function, just as in the case of the pressure above, with the result 

°! = J ** "UnM; d n = j^f^, (66) 
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where v av is the average collision frequency as determined by the loss part of the right hand 
side of the hard sphere BBGKY hierarchy [13] for the HCS 

A-K d/2 a* d ~ l f 

Vav = r J dg* 12 e(q*i2 ■ 9I2W12 ■ 9*i2 F h {9*12, x) ■ (67) 

This is a purely collisional quantity, arising from the boundary condition associated with 
hard sphere dynamics at the point of contact. 



E. Bulk Viscosity 

The bulk viscosity has forms similar to those reported above for the shear viscosity, 
although it represents resistance to volume dilation rather than shear, 

k* = lim Vl K H (s) = + lim / ds' £l G (s') . (68) 

Jo 

The correlation functions in the above expressions are shown in Appendix F to be given by 
Q H ( s ) = -(yd 2 )" 1 J dT* tr H* f (r*)e~*( C + ^M* K (T*), (69) 



tt K G (s) = -(V*^)' 1 J dT* tr H* f (T*)e V +2 /T*(r*), (70) 
and = tt H (0). Here H*^ is the projection of the integrated momentum flux, whose trace 



is 



tr H* / (r*) = tr H*(r*) - ^|j^jV - ^E* - ^iv) p* h , (71) 

where is H*(T*) is given by Eq. (43). The additional terms on the right hand side of the 
above equation, proportional to N and E*, result from projection of tr H* orthogonal to the 
eigenfunctions defined in Eq. (34). Further details can be seen in Appendix F and ref. [3]. 
The moment M* K is the same as that for the Euler transport coefficient (* u , i.e., 

N 8 

M:=Miv = -Y,<£-Q^pH(n- (72) 

r=l r 

Then, also for the associated Green-Kubo conjugate flux T*, 

t; = tt> = - (r + 1) M* K . (73) 

The instantaneous part of the bulk viscosity, Qfi, is simply related to that for the shear 
viscosity by 

n* - d -^w - ( 1 + a ^\ . f74 ) 

d ~ Ad 2 
with the average collision frequency, v av , given by Eq. (67). 
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F. Thermal Conductivity 



The expressions for the dimensionless thermal conductivity A* of the granular fluid of 
hard spheres or disks obtained in Appendix F are 

A* = lim Vl^ (s) = fio + lim / ds ' ( s ') , ( 75 ) 

Jo 

where 



n x H (s) = -(V*d)~ 1 J dT*S* f (T*)-e \ C VA4*(r*), (76) 

ft* (s) = -(Vd)- 1 j dV* S* f {T*) ■ e~"( £ ~^) S T*(r*), (77) 

and £7q = (0). In the above expressions, the projected energy flux S*f is 

S *f(r*) = S*(T*) - + plj P*, (78) 

with 5*(r*) = s*(r*;0) being the volume integrated heat flux obtained from Eqs. (D15) 
and (D29), 

r=l - r=l s^r 

(%s ■ 9rs) 2 (firs ■ G rs) firs' ( 79 ) 



As with the previous transport coefficients, the additional term proportional to the total 
momentum P* in Eq. (78) results from projection of S* orthogonal to the eigenfunctions 
defined by Eq. (34). The moment Ad* x is 

1 N 8 

= -_£,,;_. [v*Mn\ (so) 

r=l r 

and the expression of the associated Green-Kubo conjugate flux is 

r*(r*) = - (c* - f ) M* x . (81) 

The relevant hydrodynamic excitation that is subtracted from the generator of the dynamics 
is now that for A = Co/2. 

The instantaneous part of the thermal conductivity in the Green-Kubo representation 
can be simplified to 

^ = il+ iQd d+l j d °* j d 9l2 l dG 12 © (** • til) (** ■ 9l2) 

[8 (B ■ G\ 2 f + {* ■ g{ 2 f] ff ] {v\ vl vl) , (82) 
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with a* = (7*<r*. Again, the instantaneous contribution to the transport coefficient is purely 
collisional, reflecting its origin as the discontinuity in time for colliding configurations at 
contact. Also note that this is the first transport coefficient where the center of mass 
momentum of the pair, in addition to the relative velocity, occurs in the final average. 

The corresponding normal fluid expression for this transport coefficient can be analyzed 
along the lines done for the viscosities above. The Green-Kubo integrand for the time 
integral is a flux-flux correlation function as in Eq. (65), with both the direct flux S* and 
an conjugate flux, the volume integrated time reversed energy flux S*~ . This and the non- 
vanishing instantaneous part is a peculiarity of the singular hard sphere interactions and 
neither effect occurs for nonsingular forces in the elastic limit. For the granular fluid, the 
conjugate flux is further modified by its generation from the HCS. 



G. The ^-Coefficient 

The /i coefficient is a measure of the contribution to heat transport in the fluid due to 
spatial gradients in the density field. In granular fluids, density and temperature are inti- 
mately coupled through the cooling rate. At constant cooling rate, a local density variation 
gives rise to a variation in the amount of energy lost locally and hence can give rise to a heat 
flux. Since this transport coefficient is directly related to the inelasticity of the collisions, 
it is unique to granular fluids and has no analogue for normal fluids. The simplest form for 
this transport coefficient is obtained when considering the linear combination 

7 r = ^_ 2 |^A*, (83) 
am n h 

where A* is the dimensionless thermal conductivity discussed in the previous subsection. The 
intermediate Helfand and Green-Kubo representations for /J* are identified in Appendix F 

as 

J? = lim tf H (s) = Q(f + lim I ds'tt^ (s') . (84) 

Jo 

The expression of the Helfand correlation function is 

if H (s) = ~{V*d)- 1 J dT* S* f (T*) ■ e-^M^T*), (85) 
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where S* f is the subtracted heat flux given in Eq. (78). The moment is 
= [lv Q (T)] Nd J drkr 



Spi h d In Co „ 5pi h 
n- — ^— — 2— 1 



= [£v (T)] Nd n J drk* • r* 



Sn(r) d\nn ST (r) 

5pih[T\{y a }] 

Sn(r) , t 



{y a }={n,T,0} 



(86) 



{y a }={n,T,0} 

As indicated, in the last expression the functional derivative with respect to the density field 
is to be carried out at constant cooling rate. This suggests that p* is characterizing density 
gradients under conditions of constant cooling rate. 
The Green-Kubo correlation functions are 

n£ = % (o) , (87) 

fig ( s ) = -(V^dy 1 J dT* S* f (T*) ■ e- sZ Y^(r*), (88) 
with the associated flux T^* having the expression 

= -fM* (89) 

The hydrodynamic mode subtracted from the generator C* in this case has the eigenvalue 
A = 0. Unlike the previous transport coefficients discussed so far, the moment and the 
flux Ti- are still given implicitly in terms of the local HCS state. Further analysis requires 
determination of the density dependence of the inhomogeneous HCS, i.e. that in an external 
conservative field. This is beyond the discussion here and remains an open problem for 
applications of these representations for the p* coefficient. 

It is shown in Appendix G that this transport coefficient vanishes in the elastic limit, as 
expected. However, that result is specific to the choice of an equilibrium Gibbs state for 
the reference ensemble. More generally, for non-equilibrium reference states, it is p* ^ 0. 
As noted above, the linear combination p* = p* — 2 1 A* is the transport coefficient 
measuring the effects of density gradients on the heat flux when the cooling rate instead of 
the temperature is held fixed. The Helfand moment Ai^- is the space moment associated 
with ^;(r*;0) - 2f^^(r*;0). It is shown in Appendix C that this same combination 
determines the eigenfunction of the specific hydrodynamic mode solution of Eq. (34) with 
eigenvalue A = 0. Hence, it appears that in some respects the pair S(o, 5n are more natural 
independent variables than 5T, 5n. 
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IV. EVALUATION OF TRANSPORT COEFFICIENTS 

The (intermediate) Helfand and Green-Kubo expressions for the various transport coeffi- 
cients obtained above, are formal results in that the full N particle problem has not yet been 
solved. Instead, the linear response method gives a valuable exact and direct relationship 
of the macroscopic properties to the microscopic statistics and dynamics. Their utility is 
determined by the fact that they are the ideal setting to explore controlled analytic approx- 
imations and exact numerical evaluations. The transport coefficients have been obtained 
in the form of stationary time correlation functions. This suggests the possible general- 
ization to granular fluids of the extensively developed many-body tools for the analysis of 
time correlation functions of a normal fluid, such as short time expansions, mode coupling 
theories, and formal kinetic theories. Also, for normal fluids such expressions have proved 
particularly suitable for evaluation by MD simulation to extend the studies of hydrodynamic 
descriptions to otherwise inaccessible domains in the density of the fluid. In this section, 
this potential for future studies is illustrated in two directions. First, the applicability of a 
kinetic theory of time correlation functions is outlined and discussed. Then, a scheme for 
implementing the stationary representation of the dynamics considered in this work in an 
MD simulation is discussed. 

A. Kinetic Theory 

The results for the various transport coefficients obtained in the present work are ex- 
pressed in terms of time correlation functions over an N particle distribution. However, an 
equivalent expression involving two particle time-dependent reduced functions is possible as 
well. This reduced representation serves as the appropriate starting point for construction 
of a kinetic theory for the correlation functions, and application of formal cluster expansion 
techniques. Consider the Helfand representation for a generic transport coefficient to* as 
given by Eqs. (30) and (32) 




(90) 
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The projected source or flux, F* s 'f (T*), is the sum of single particle or pair functions in all 
cases, so it has the form 

N N N 

F *sj (n = Fl (x* r ) + £ x> 2 w> *:) > ( 91 ) 

r=l r=l Sy^r 

with x* = {qr*, v*} denoting the one particle phase point for particle r, and F\ and F 2 being 
one and two particle functions of the phase point, respectively. Define a set of "reduced 
correlation functions", M*^ (x*, ...,x* m , s) by 

(xl,...,x* m ,s)= ™ J dx* m+1 ...J dx* N e-< r - x )M*(r*). (92) 

Then the generic expression for the transport coefficient (90) becomes 

zu* = \imV*~ 1 J dx* 1 F 1 (x* 1 )M* {1) (xl,s) + J dx\ J dx* 2 F 2 (x*, x* 2 ) M* {2) {x\, x* 2 , s) . 

(93) 

In this way, the full N body correlation function is mapped onto one that involves only the 
one and two particle reduced correlation functions. 

It follows by direct differentiation with respect to s of Eq. (92) that the functions Ai*^ 
obey a hierarchy of equations, 

= E / dx *m + i r (r, m+1) Af ( ™ +1) {x{, x* m+1 , s) , (94) 

r=l •* 

where C is the same operator as C but now defined for only m particles, and the T (r, s) 
operator is the dimensionless form of the operator defined in Eq.(A24). This hierarchy is 
formally the same as the dimensionless BBGKY hierarchy of the reduced distribution func- 
tions representing the non-equilibrium state of the fluid [13, 15]. Only the initial conditions 
are different. 

Formally, a kinetic theory represents a closure of the hierarchy at the level of the first equa- 
tion, through the identification of a representation for the two body function Ai*^ (x*, x 2 , s) 
as a functional of the one body function Ai*^ (x*, s). With such a functional relationship, 
the first hierarchy equation becomes a closed, autonomous equation for M*^ (x*, s). This 
equation is called a kinetic equation, and the functional relationship defines the kinetic the- 
ory. It is important to recognize that this concept of a kinetic theory does not imply any a 
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priori restrictions to low density, weak inelasticity, or absence of correlations. Such restric- 
tions occur only in the justification of specific choices for approximate closure functionals. 

To obtain a functional relationship for M.*^ (x\,x* 2 , s), it is postulated that there exists 
a functional expressing the two body reduced distribution function associated with the local 
HCS in terms of the one body distribution function, 



fih (x U X 2 ,t) = \x U X 2 ,t\f£ 

where the reduced local HCS distributions are defined by 



(2) 



? (1) 



Am) 
Jlh 



X\ . . . , X m , t) — 



dx 



m+l 



dx N e tL pi h (T). 



(95) 



(96) 



(n — m)\ 

The formal construction of a functional like in Eq. (95) for normal fluids, where pih is replaced 
by the local equilibrium ensemble, is a well-studied problem, e.g. inversion of formal cluster 
expansions for the reduced distribution functions. This is the first point at which extensions 
of kinetic theory to granular fluids can be explored. From the definition of Ai* in Eq. (33), 
it follows that Eq. (95) implies a similar but linear functional relationship for M*^ , 



with 



K{x\,x* 2 ,x*', s) = [lv (T h ) 



(i) 



x\,x 2 ,t\fi h 



8fS? (^t) 



(97) 



(98) 



{y a }={n h ,T h ,0} 

The linear kernel K is a functional derivative of the two body distribution function with 
respect to the one body distribution, evaluated in the homogeneous limit. This latter eval- 
uation leads to the linearity of this functional. The pre-factor in the above expression just 
transforms to dimensionless variables as in Eq. (19). Substitution of this form for the two 
body correlation function in the first equation of the hierarchy leads the kinetic equation 

(d s - A + - 1 (s)) (x{, s) = 0, (99) 

with the formal linear "collision operator" 

J (s) M* (1) (xl, s) = J dx* 2 T* (1, 2) J dx*'K (xl, x* 2 , x*', s) M* {1) (x*', s) . (100) 

This is the kinetic equation governing the dynamics of M*^. Since this quantity also deter- 
mines M*^ through the functional relationship given by Eq. (97), the transport coefficient 
can be calculated from Eq. (93). 
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The above is exact but formal, and the aim here is only to provide a flavor for the potential 
for developing a practical kinetic theory for time correlation functions. Details of this method 
and the explicit evaluation of the various transport coefficients an Enskog-like approximation 
will be given elsewhere [16]. The latter is obtained by estimating [xi, X2, t\f^~\ ~ 
[xi,X2, 0\f^] , which is exact at t — 0, and furthermore neglecting velocity correlations. 
The resulting kinetic theory for the correlation functions leads to cooling rate, pressure, and 
transport coefficients identical to those from the Chapman-Enskog solution to the Enskog 
kinetic theory [17, 18]. At low density, these results agree with those from the Boltzmann 
equation [19]. 



B. Molecular Dynamics Simulation 

The formal expressions for the various transport coefficients derived in this work, are 
given in the dimensionless representation for the phase space dynamics that leave the HCS 
ensemble stationary and defined by the thermal velocity associated to the HCS. But for 
practical purposes, this is a complicated representation as it must be solved self consistently 
with the cooling equation. It is then convenient to consider a different scaling, in which the 
function Vo(t) is replaced by a known function chosen in a appropriate way. A useful choice 
is defined by 

n* ft 

q r =q r = — , v r =T^v r (101) 

where £0 is an arbitrary time-independent dimesnionless frequency. The rational for his 
choice and some more details are given in Appendix H. The particle dynamics in the new 
variables consists of an acceleration streaming between collisions, 

dq**(r) 



dr 



v**(t), (102) 



^ = |«-M, (103) 

while the effect of a collision between two particles r and s is to instantaneously alter their 
relative velocity according to the same rule as in the original scale. This invariance of the 
interation law is another consequence of the singular character of hard particle collisions. 

The quantity £0 an d the dimensionless cooling rate Q are related through the steady state 
temperature obtained in the simulations in the scaled variables [20, 27]. The relationship is 
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given in Eq. (G15). This method has already been implemented in event-driven molecular 
dynamics simulation to measure the self-diffusion transport coefficient of a granular fluid 
[20]. That is a simpler case than those considered here, since the correlation functions 
appearing in its Helfand and Green-Kubo expressions are not moments of the conjugate 
densities, but simply the velocities of the particles. 

The low density limits of the Green-kubo expressions for a granular gas have been eval- 
uated by means of the direct simulation Monte Carlo method in refs. [21] and [22], also by 
using the steady representation introduced above. The conjugate fluxes were determined 
from the simulation themselves. Although the technical problems are much harder beyond 
the low density limit, this shows the way to evaluate the general expression derived here at 
arbitrary densities. 

V. DISCUSSION 

The parameters of Navier-Stokes hydrodynamics for a system of smooth, inelastic hard 
spheres with constant restitution coefficient have been given formally exact representations 
starting from the underlying non-equilibrium statistical mechanics. General linear response 
functions for small initial spatial perturbations of the HCS are given by Eq. (15). These are 
stationary time correlation functions, composed of densities of the hydrodynamic fields and 
corresponding conjugate densities defined by functional derivatives of the local HCS with 
respect to these fields. They are similar to the response functions for a normal fluid, where 
the local HCS is replaced by a local equilibrium ensemble. Such response functions provide 
the basis for exploration of a wide range of phenomena in granular gases at the fundamental 
level of statistical mechanics, beyond the specific case of transport coefficients considered 
here [23]. 

In the long wavelength, long time limit these response functions exhibit the hydrodynamic 
behavior of the system. This is the analogue of the Onsager regression property for normal 
fluids [24, 25] . This limit has been extracted by defining a general transport matrix for the 
response functions, and expanding it to order k 2 to identify the Navier-Stokes parameters in 
terms of corresponding HCS averages or stationary time correlation functions. The cooling 
rate is given as the HCS average rate of change of the global energy, and the pressure 
is expressed as the average fluctuation in the global momentum flux. Such averages can 
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be studied directly by MD simulation. Furthermore, they have been expressed in terms 
of the reduced two particle distribution function for pairs of particles at contact. This 
representation is suitable for inquiry about new features of granular fluids, such as two 
particle velocity correlations. 

The transport coefficients obtained in this way have been considered in two represen- 
tations, the (intermediate) Helfand and Green-Kubo forms. The former is the long time 
limit of a correlation function constructed from a global flux associated with the energy 
or momentum densities, and a space moment of one of the conjugate densities. This is 
quite similar to the corresponding Helfand results for a normal fluid, if the HCS ensemble 
is replaced by the equilibrium Gibbs ensemble. There is an important additional difference, 
however, in the generator for the time dependence. The dimensionless generator for hard 
sphere trajectories L* is replaced by £* of Eq. (22), which contains an additional velocity 
scaling operator. The latter represents the cooling of the reference HCS and it is essential 
for the stationary representation given here. In addition, the time correlation functions for 
transport coefficients have the further modification £ — > £ — A, where A is one of the 
eigenvalues of the hydrodynamic transport matrix JC* hyd (0). These eigenvalues describe the 
dynamics of homogeneous perturbations to the HCS. Since the transport coefficients mea- 
sure only spatial correlations, it is perhaps not surprising that this homogeneous dynamics 
does not occur in their representation. These same eigenvalues occur in the exact spectrum 
of £* , so this generator £* — ]C* hyd (0) has invariants (zeros of the spectrum). In the elastic 
limit, these invariants are proportional to the usual global conserved total number, energy, 
and momentum. 

The equivalent Green-Kubo expressions have a time independent part plus the familiar 
time integral of a flux-flux correlation function. The time independent part is absent for 
normal fluids with conservative, nonsingular forces. More generally, for either singular forces 
(e.g., hard spheres) and/or nonconservative forces, this contribution is nonzero. For normal 
fluids, it can be evaluated exactly and represents the high density contribution to the trans- 
port coefficients in the Enskog kinetic theory. Here, it can be reduced to an average over 
the two particle distribution function, just as described for the cooling rate and pressure. 
The flux-flux time correlation functions in the Green-Kubo integrands are constructed from 
a global flux of the energy or momentum densities and a global flux associated with the 
conjugate densities. The dynamics is generated by £* — JC* hyd (0) as described above, so the 
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associated invariants must not occur in order for the time integral to exist in the long time 
limit. In fact, the expressions obtained have the fluxes orthogonal to these invariants. This 
is the analogue of "subtracted fluxes" for a normal fluid [13]. 

There are additional transport coefficients associated with spatial gradients of the hydro- 
dynamic fields contribution to the cooling rate. These also have Helfand and Green-Kubo 
forms as described above, but with the global flux of energy replaced by the volume in- 
tegrated source in the microscopic balance equation for energy. One of these transport 
coefficients, ( u , appears at the Euler level hydrodynamics and gives a correction to the 
pressure in the coefficient of V ■ U. The other two give corrections to the macroscopic heat 
flux proportional to VT and Vra. These vanish for a normal fluid. 

As with normal fluids, the value of such formal expressions is direct access to the trans- 
port coefficients before approximations are introduced. Often this avoids (or postpones) 
conceptual difficulties associated with the particular many-body method used for evalua- 
tion. Here, the two most common approaches for normal fluids have been briefly discussed 
for extension to granular fluids. The first is a kinetic theory for time correlation functions. 
It is inherently simpler than constructing a corresponding kinetic theory for the reduced 
distribution functions, since the latter are nonlinear while that for the correlation functions 
is linear. Furthermore, the special solution required for hydrodynamics (normal solution 
and Chapman- Enskog construction) has already been incorporated implicitly in the linear 
response analysis here. The necessary functional assumption for a kinetic theory of gran- 
ular time correlation functions has been identified as a property of the local HCS. In this 
restricted context, it should be possible to explore the validity or failure of such an as- 
sumption, although the detailed investigations for normal fluids need to be reconsidered. It 
is shown elsewhere [16] that an ad hoc neglect of velocity correlations, something that is 
justified in the case of normal fluids, leads to transport coefficients identical to those from 
the Enskog kinetic theory. The formalism developed in this paper and the preceding one, 
provides the basis for critiquing that theory and correcting it for granular fluids. 

The second approach for evaluating the formal expressions given here is MD simulation. 
In this case, the extension of approach from normal fluids is not direct either. The conjugate 
moments (Helfand representation) or the conjugate fluxes (Green-Kubo representation) are 
not known explicitly, but rather given by derivatives of the HCS ensemble. Thus the simu- 
lation of these expressions requires new method for developing properties of the HCS from 
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the simulation itself. 
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APPENDIX A: HARD SPHERE DYNAMICS 

In this Appendix, the details of the dynamics and statistical mechanics of a system 
of hard particles is given. Consider a system of N mono-disperse smooth hard spheres 
(d = 3) or disks (d = 2) of mass m and diameter a. A complete specification of the 
initial state of the system involves knowing a point in the 2iV<i-dimensional phase space 
T = {q r , v r ;r = 1, . . . , N}, that gives the position q r and velocity v r of each particle r of 
the system. At a later time t, the state will be characterized by another point in phase space, 
r t = {q r (t), v r (t); r — 1, . . . , N}. The dynamics of the particles consists of free streaming 
(straight line motion along the direction of the velocity), until a pair of particles, r and s, are 
at contact, at which time their velocities v r ,v s change instantaneously to v' r ,v' s according 
to the collision rule 

' _ 1 + Q y - \~ 

V r V r \Qrs ' 9rs)9rsi 

1 + a 

v' s = v s + -^—(q rs ■ g rs )qrs, (Al) 

where q rs is a unit vector along the relative position q rs = q r — q s at contact and a is the 
coefficient of normal restitution, defined in the interval < a < 1 and that will be taken 
here as constant, i.e. velocity independent. In terms of the relative velocity g rs = v r — v s 
and the center of mass velocity G rs = (v r + v s ) /2, the above collision rule reads 

G' rs = G rs , g rs = 9rs ~ (1 + a) (q rs ■ 9rs) q rs - (A2) 
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The energy loss on collision is 



(A3) 



It is seen that a — 1 corresponds to elastic collisions. Subsequent to the change in relative 
velocity for the pair {r, s}, the free streaming of all particles continues until another pair is at 
contact, and the corresponding instantaneous change in their relative velocities is performed. 
The sequence of free streaming and binary collisions determines a unique trajectory in phase 
space, T t , for given initial conditions. The collision rule is invertible so the trajectory can 
be reversed, although with the inverted collision rule ("restituting" collisions). 

The statistical mechanics for this system [26] is comprised of the dynamics just described, 
an initial macrostate specified in terms of a probability density p(T), and a set of microscopic 
observables (measurables) . The macroscopic variables of interest are the expectation value 
for the observables at time t > for an initial state p(T). Given an observable A(T), its 
expectation value is defined by the two equivalent forms 



In the above expression, L is the generator of the dynamics of phase functions and L is the 
generator of dynamics for distribution functions. The explicit forms of these generators will 
be derived in the following. 

1. Generator of Trajectories 

The hard sphere dynamics described above is characterized by piecewise constant veloc- 
ities that change instantaneously (and discontinuously) at the time of collision. This allows 
the generator of trajectories to be derived using geometric arguments [10]. Consider first a 
system of two inelastic hard spheres or disks. The particles move freely until they eventually 
contact, at which time their velocities change instantaneously according to Eq. (Al) above. 
Suppose that the two particles, named 1 and 2, respectively, collide at time r. This time is 
determined as a function of the initial separation q 12 and relative velocity gi 2 from 



with a being an arbitrary vector of modulus a. Of course, if this equation does not have a real 
and positive solution for r, the particles do not collide. Then, for < t < r the trajectory is 




(A4) 



(A5) 
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given by free streaming with the initial velocities, i.e. T (t) = T (t) = {q t + v r t, v r ; r — 1, 2}. 
For t > r, it is T (t) = T'(t) = {q r + v r r + v' r (t — r) , t^,; r = 1, 2}, where the velocities u£ 
are determined by the collision rule (Al). Therefore, the value of any phase function A (T) 
at time t can be given compactly expressed as 

a [r(t)\ = {i - e [t - r (r)]} a [To (*)] +e[t-r (r)] a [r' (*)] , (A6) 

where Q(x) is the Heaviside step function, defined by Q(x) = for x < and 6(x) = 1 for 
rr > 0. Differentiation with respect to time of the above equation leads to 

= ± Vr (t) . «*ra + s lt _ T (r)1 M r wl - , |ro(( )]} 

2 

= E^(t)-^|^ + 5[t-r(r)](6 12 -l)A[r(t)]. (A7) 

Here, &i 2 is a substitution operator that changes the velocities into their postcollisional values 
as given by Eq. (Al), 

b 12 X v 2 ) = X (6i2«i, &12V2) = X K, v 2 ) . (A8) 

Of course, the vector q\ 2 to be used in the collision rule is <7i2(t)/<?i2(t)- 

It is convenient to substitute the delta function involving the time t in Eq. (A7) by a more 
geometrical quantity. To begin with, notice that not every solution of Eq. (A5) corresponds 
to a physical collision. In fact, it is easily realized that for every contact representing a 
physical collision, there is another solution to Eq. (A5) corresponding to the time at which 
the two particles would separate if they were allowed to stream through each other. The 
physical value of r can be identified as that for which 912 (t) g\ 2 < 0, while for the unphysical 
solution of Eq. (A5) it is <Zi 2 (t) • 912 > 0. Next, take into account that the solutions of Eq. 
(A5) are the same as those of qi 2 (r) — a — 0. Hence 

dg 12 (t) 
dt 

© [-912(0 • 012(f)] 5 [q 12 (t) - a] \q 12 (t) ■ g 12 {t)\ . (A9) 



5[t-r(T)] = e[-q 12 (t)-g 12 (t)]5[q 12 (t)-a] 



Use of this in Eq. (A7) yields 

dA[T(t)] A , , dA[T(t)} rr . , ^ . 

ft = E ^ (0 • J (V + * ^ (0 - 6 [~gl2 (t) ■ gi2 (*)] \qi2 (t) • 012 (01 

x(fe 12 -l)A[r(0]. (A10) 
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The generator of trajectories is defined by 

A [r (t)] = e Lt A (r) , (All) 

so L is identified as 

2 r, 

L = E^-7^T + T ( 1 ' 2 )' ( A12 ) 

r=l 9r 

with the binary collision operator T(l, 2) given by 

T(l, 2) = 5 (g 12 - a) (-q 12 • g 12 ) \q 12 ■ g 12 \ (b 12 - 1) . (A13) 

This expression can be rewritten in a form that is often more suitable for explicit calculations, 
by employing the relation 



da S(q 12 -a) = a-^5(q 12 - a) , (A14) 

where da is the solid angle element associated to a. In this way, it is obtained, 

T(l, 2) = a d - 1 J da S(q 12 - <r)Q(-a ■ g 12 ) \a ■ g 12 \ (b 12 - 1). (A15) 

The generator of trajectories for the system of N particles considered here, entails the 
additional assumption that only binary collisions occur between particles, i.e., there are no 
collisions involving three or more particles simultaneously at contact. Then, for a system of 
N particles, 

N „ N N 

L = H V r- W + - 2 Y.Y. T ^ S )- (A16) 

r=l r sjtr 

2. Adjoint Dynamics 

Next, consider the generator L defined by the second equality in Eq. (A4). It can be 
identified from 

J dTA(T)Lp[T] = - J dT [LA(T)\ p(T) 

N . „ N N . 

= -E/ rfr ^ r )^-^-A(r)--^^ / dTp(T)5( qrs -a) 

r=l r s^r 

X0 (-q rs ■ g rs ) \q rs ■ g rs \ (b rs - 1) A (T) . (A17) 
Define the inverse, ft" 1 , of b rs by b~}b rs = b rs b~} = 1. Equation (A2) gives directly 

Ks9rs = 9rs = 9rs (firs ' 9rs) <7rs ( A18 ) 
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A useful identity is given by 

J dvx (r) b rs Y (r) = J d ip-}v) x ip-}v) y (r) = cr x J dvx (b^r) y (r) , (Ai9) 

for arbitrary functions X(T) and Y(T). The first equality is obtained by changing integration 
variables from (v r ,v s ) to (b ). The factor a 1 is the Jacobian for this change of 

variables. Also, 

Ks (Qrs ■ 9rs) = -U~ l q>rs ' 9rs (A20) 

Equation (A17) therefore can be rewritten as 



N „ „ AT 



y rfrA(r)Lp(r) = J^J drA(r)v r ■ ^-p(r) y ds r -y dv r y rfr (r) v rP (r) a (r) 

_ 2 / rfr ^ ( r ) 5 (<?rs - O") Wrs ■ 9rs\ 

x [0 (£ rs • g rs ) a-%, 1 - (-q„ • g rs )] p (T) . (A21) 

An integration by parts has been performed in the first term, leaving a surface integral 
for which boundary conditions must be specified. Then dS r is the surface element vector 
associated to q r and dT^ is the element of volume of the 2(N — l)d phase space obtained by 
eliminating q r and v r in the original one. This surface term identically vanishes when, for 
instance, periodic boundary conditions are considered for the system and A(T) is a compact 
function of the positions. In the following it will be always considered that this is the case. 
The adjoint Liouville operator is formally identified from the other non-surface terms as 

N n 1 AT N 

i=r r s^r 

where the new binary collision operator is 

T(r, s) = S(q rs - a)\q rs ■ g rs \ [0 (q rs ■ g rs ) a~ 2 b~} - (-q rs • grs)] • (A23) 
The integral representation of this operator, similar to Eq. (A15), is 

T(r, s) = G d - X y d$e(a ■ g rs )\a ■ g rs \ [a- 2 5(q rs - (r)b;} - 5(q rs + a)] . (A24) 
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3. Time Correlation Functions and Time Reversed Generator 



The response functions in statistical mechanics take the form of time correlation func- 
tions over the macrostate under consideration, and can be written in two equivalent ways, 
corresponding to each of the two representations in Eq. (A4), 

C AB (t) = JdT [e tL A(T)} B(T)p(T) = J dT A(T)e' tT [p(T)B(T)] , (A25) 

where A(T) and B(T) are two phase functions. A third representation can be identified in 
the form 

C AB (t) = JdT A(T) [e- tL -B (T)} e~ a p (T) . (A26) 
Comparison of Eqs. (A25) and (A26) shows that the operator L_ must verify the relation 

L(pB) = (Lp)B + pL_B. (A27) 

It follows that 

N BE 1 N N 
pL_B = p vr ■ g- - 2 E E s ) (P B ) ~ BT ^ S M ■ (A28) 

r=l r=l Sy^r 

Next, consider the action of T(r, s) on pB in detail using the definition (A23), 

T(r, s) (pB) = 5(q rs - a)\q rs ■ g rs \ [Q (q rs ■ g rs ) a' 2 (6"V) ip~}B) 
-0 (-q rs ■ g rs ) pB] 
= BT(r,s)p 

+S(q rs - a)\q rs ■ g rs \Q (q rs ■ g rs ) (a~ 2 6~V) (b^ 1 - l) B 
= BT(r, s)p + p5{q rs - a)\q rs ■ g rs \Q {q rs ■ g rs ) (6" 1 - l) B. (A29) 

In the last equality, use has been made of the identity 

5(q rs - a)\q rs ■ g rs \0 (q rs ■ g rs ) a~ 2 b~^p = 5(q rs - a)\q rs ■ g rs \0 (q rs ■ g rs ) p, (A30) 

valid for hard sphere or disk distribution functions [27]. This identity follows from the 
requirement that the flux for a pair of particles at contact with relative velocity g rs on the 
pre-collision hemisphere, should be the same as that for g" s on the post-collision hemisphere. 
See ref. [27] for further details and applications. 
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Use of Eq. (A29) into Eq. (A28) gives 

dB N N 

pL_B = p^V r ■ | ^2^2 6 ((Irs - vMrs ■ 9rs\&(qrs ' 9rs)(Ks ~ X )- ( A31 ) 

r r=l Sy^r 



This in turn gives the identification 

N r, N N 

L - = J2 V r-g--^J2 T -^^ ( A32 ) 
r=l r=l s^r 

with a third binary collision operator TL(r, s) defined as 

T_(r, s) = <5(g rs - cr)e(g rs • <? rs )|q rs • g„)| - l) 

= a 4 ' 1 J da Q(a-g rs )\q rs -q rs )\5(q rs - (7)^-1). (A33) 

It is a simple task using standard techniques to verify that Eq. (A27) implies 

e - t7: (pB) = (e~'V) e- tL -B, (A34) 

and Eq. (A26) follows. For the particular case p = ph, the HCS ensemble, it is 

e- {L Ph [T;T(0)]=p h [T;T(t)} (A35) 

and, therefore, the correlation functions defined in Eq. (A26) can be expressed in the form 

C AB {t) = J dY A(T) [e- tL - B(T)] p h [T; T(t)} . (A36) 

It is seen that L_ is the same as the trajectory generator identified in Eq. (A16), but 
with the restituting collision rule in place of the direct collision rule. Therefore, L_ is the 
generator of time reversed trajectories in phase space. 

APPENDIX B: DIMENSIONLESS REPRESENTATION 

The time dependence of the reference HCS ensemble has the scaling form indicated in 
Eq. (4). This means that the time dependence can be removed entirely by a change of 
variables, leaving a universal dimensionless function of the scaled velocities. It is useful to 
introduce these same variables more generally for other states as well, to partially account 
for collisional cooling. The advantages will be apparent in the final results. The set of 
dimensionless variables are chosen to be 

q r = -, Vr = W y ds=—dt, (Bl) 
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where v (t) = v [T h (t)] is defined in terms of the temperature of a HCS reference state 
having the same initial total energy as the actual system under consideration. The Liouville 
equation then becomes 

() +rV(r, a ) = o, ( fi 2) 



ds 

with the reduced distribution function p*(T*, s) defined by 

p*(T*,s) = [£v (t)] Nd p(T,t), (B3) 
being T* = {qr*, t>*; r = 1, . . . , N} and the operator C* given by 



N 

2 ^-f dv, 

r=l ' 

Here, 



£V = ^V + f £^-KV*)- (B4) 



L* = -^ L= [ L ]r=r* (B5) 
and Q is the dimensionless cooling rate 

Co " ■ (B6) 

— * — * 

The first term, L , of the generator C for dynamics in this dimensionless Liouville equation, 
generates the usual hard sphere trajectories in scaled variables. The second term rescales 
the velocities along the trajectories, to represent their change due to the average cooling 
associated with the HCS. A first advantage of this representation is the fact that the HCS 
ensemble is a stationary solution of the Liouville equation, 

Tp* h = 0. (B7) 

Consider a phase function A(T) whose dimensions scale out as a factor with the above 
change of variables, i.e., 

A(r) = c A [v (t)]A*(F), (B8) 

where [f (t)] contains the dimensions of A(T) so that A* (T*) is dimensionless. The 
ensemble average of A (T) at time t is 

(A (t)) = J ' dTp (r, t) A (r) = c A [v (t)} (A* (s))* . (B9) 

The dimensionless average is then given by 

(A* ( s )>* = ^ = [ dr* \e~ sT P * (r, 0)1 a* (r*) . (bio) 

ca [vo (t)\ J L J 
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Adjoint operators can be introduced just as in Appendix A for equivalent representations. 
For example, 

(A*(s))* = J dT* p* (T*,0)e sC *A* (T*) (Bll) 



with 



r=l 



L * = ilt) L = [i|r - ' (B13) 

To illustrate the utility of this formulation, consider an average over the HCS. Then Eq. 
(BIO) gives (A* (s))* = (A* (0))*, since is stationary. Using Eq. (Bll), the time derivative 
at s = is seen to be given by 

/W)^.(r-)=0. (B14) 

Suppose now that A* be an arbitrary differentiate function of the scaled total momentum 
of the system P* = ^2 r v*, i.e. A* = A*(P*). By momentum conservation, L*A*(P*) = 
and Eq. (B14) becomes 

| J dT* pi (P) < ■ ( p l = f / dr* p* (P) P* • ^A* (P*) = 0. (B15) 



Since this holds for any function A*(P*), it follows that pjj(r*) must be proportional to a 
delta function in the total momentum, 

P i(n = 5(p*)pi(n- (Bi6) 

This is the result quoted in Eq. (5) in the text. 

The dimensionless form for correlation functions defined in Eq. (A25) are obtained by 
using the representation in Eq. (Bll), 

^" % A Jh M r / r [e-A-(r.)] B - ( r>. r ,o). (BIT) 

The adjoint representation follows directly as well, 

C* AB (s) = J dT* A* (P) e~ sT [p* (P, 0) B* (P)] . (B18) 
Finally, the representation in Eq. (A26) gives the third equivalent form: 

C* AB (s) = J dT* A*(T*) [ e - sZ V (r*, 0)] e- sC -B*(T*), (B19) 
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with 

t* N r) 

C_ = L*_ + % V v* ■ (B20) 
2 ^-f dv* 

and 

L*_ = = [L_] r=r , . (B21) 

For the particular case of the HCS, in Eqs. (B17)-(B19) it is 

e~ sT pi (r* , 0) = pl{Y\ 0) = pj(r), (B22) 

due to the stationarity of the distribution of the HCS in the reduced scales. 

The results of this Appendix, starting from hard sphere dynamics, agree with those 
obtained in ref. [3], starting from more general interaction potentials and taking the scaling 
limit (see Sec. 8 in ref. [3]). 



APPENDIX C: HYDRODYNAMIC MODES 



In this Appendix, the eigenvalues of the hydrodynamic transport matrix JC* hyd (k*), whose 
expression is given by Eqs. (11)- (13), are presented and some comments made on their 
implication for transport in a granular fluid. Similar results, but restricted to the low density 
limit have been reported in [28]. The eigenvalues of the matrix ]C* hyd are the solutions j a 
to the cubic equation 
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k* 2 = 0, 



4 V^ lnn ^ Ph, 

plus the (d — l)-fold degenerated shear modes 



7± 



_(i + r] * k *2 



(Cl) 



(C2) 



If the limit a — > 1 is taken for this equation, £q — > and the solutions to order k* 2 give 
the familiar hydrodynamic modes associated with normal fluids: the two propagating sound 
modes, the heat mode and the d — 1 transverse shear modes [13]. But, when the solution to 
the above equation is considered at finite a, the modes to order k* 2 are quite different, 



Co* \dlnn h 



ln^ ) < - 

Ph 



(C3) 
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All eigenvalues are real and hence there are no propagating modes. Also, the limit a — > 1 is 
singular so these modes in this limit do not represent the familiar hydrodynamic excitations 
of a normal fluid. The drastic difference in the nature of the hydrodynamic modes obtained 
as the elastic limit of the above eigenvalues, is due to the non analyticity of the eigenvalues 
and eigenvectors about the point a = 1 and k — 0. Close to the elastic limit, both Q oc 
(1 — a 2 ) and k are small parameters, and the type of modes obtained depends on how these 
parameters approach zero [28]. This is an indication of the fact that the inelasticity, even 
when small, gives rise to drastically different transport in the fluid. For the purposes at 
hand, attention is restricted here to the q^I forms of these modes. 
There exists a critical wavelength k*f defined by 



a* 

7, *c _ SQ 



(C6) 



such that for k* < k*f the shear modes become unstable. Similarly, there exists a threshold 
wavelength for the 7y mode such that it becomes unstable as well. This means that the 
homogeneous state characterized by these hydrodynamic equations, is unstable to sufficiently 
long wavelength perturbations that excite these modes. This instability of the HCS has been 
well established in the literature [29, 30]. This implies physically that the response due to the 
unstable modes grows until such time as the linear theory breaks down, and further analysis 
of the dynamics has to be carried out using the full non-linear hydrodynamic equations. It 
is still an open question as to the nature of the final state [31]. 

Further insight into the nature of the hydrodynamic response of this fluid can be obtained 
by looking at the corresponding eigenvectors. To lowest order in k these are found to be 
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Here, = and i — 1 = 1, for d — 2, while 







W 

and i — 1,2, with 1 = 



(C7) 




and 



2 = 




for d = 3. 



The first of these modes is excited by the condition 



ST* = -2^^- 5n* (C8) 
am n/t 

and zero flow velocity. This can be interpreted as follows. The cooling rate has the form 
(o(nh,Th) = T^ 2 ( (n h ). It then follows that this condition for exciting the first mode 
corresponds to variations in the temperature and density that leave the cooling rate constant. 
The second mode in Eq.(C7) is due to a temperature perturbation at constant density and 
also zero velocity, while the third one is due to a longitudinal velocity perturbation at 
constant temperature and density. The last d — 1 modes are the response to a transverse 
velocity perturbation, again at constant temperature and density. 



APPENDIX D: MICROSCOPIC CONSERVATION LAWS 

The aim here is to derive the conservation laws (more precisely, balance equations) for 
the microscopic observables of interest, namely, the number density Af(T; r), energy density 
S(T;r), and momentum density Q(Y;r). Their ensemble averages give the hydrodynamic 
fields. One set of balance equations is associated with the dynamics for t > 0, whose 
generator is the L operator given by Eq. (A16). Another set of balance equations describes 
the dynamics for t < 0, as generated by L_ defined in Eq. (A32). Both dynamics lead to 
forms of fluxes that appear in the Green-Kubo expressions for both elastic and inelastic hard 
sphere transport coefficients (see appendix H). 

Let {A a (T;r,t)} denote the set of phase functions {N,£,Q}. Consider first t > 0. The 
phase functions obey the equation 

dA a (r;r,t) 



dt 
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L(T)A a (T;r,t) = 0, (Dl) 



where L(T) is given in Eq. (A16). The aim here is to evaluate the action of the L operator 
on the density A a to identify a balance equation of the form 

dA(T;r,t) 



dt 



+ V-j a (T;r,t) = -w a (T;r,t). 



(D2) 



Here, j a is the flux associated with the density A a , and w a is a source that signifies a local 
loss contribution in this density that cannot be expressed as a divergence. It is non zero 
only in the equation for the energy density. 
The number density is defined as 



N 



A/-(r ; r) = ]T5(r-g r ), 



r=l 



SO 



N 







N 



LMiV; r) = J2 v r ■ ~n ^ (r ~ Qr) = - V • ^ v r 5 (r - q r ) . 

Oq 



r=l 



r=l 



This gives the microscopic continuity equation, 

dtf(T;r,t) | v Q(T;r,t) _ Q 
dt m 

where the number flux density is seen to be proportional to the momentum density 

N 



0{T;r) = ^2mv r 5 (r - q r ) 



(D3) 



(D4) 



(D5) 



(D6) 



r=l 



Consider now the equation for this density. It is 

_d_ 
dq 



N q N N 

LQ(T;r) = ^ mv r v r ■ — — 5 (r — q r ) + — ^ ^ T (r,s) [v r 5 (r — q r ) + v s 5(r — q s )\ 



r=l 



-V 



r=l s^r 



N , s N N 

Er/ \ 1 + « Wlff v^v^ r/ \ r\ / ~ \ 

mv r v r 6 (r - q r ) H 2^X^ d ^ q ' rs ~ a )® • 9rs) 



r=l 



r=l s=£r 



X 



\q rs -9rs\ 2 QrsQrs / dj 5 (r - q r + ^q r 
Jo 



In the second equality, use has been made of the relation 



(b rs - l)v r = -(b rs - l)v s = 



1 + a 



(.Qrs ' Qrs) Qr 



(D7) 



(D8) 



and the identity 



f 1 d 

8( r ~ Qr) - $( r ~ Qs) = ~ dl -k-S (r -q r + -yq rs ) = -Qr 

Jo a l 



/ d-y5(r -q r + 7<? rs ) . 
Jo 

(D9) 
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Equation (D7) gives the conservation law for the momentum, 

d 



dt 



g(T;r,t) + V -h(r ; r,t) = 0, 



(D10) 



with the tensor momentum flux density h identified as 

N s V \ 

1 



h(r;r) = ^mv r w r 5(r — q r ) H (q rs - a) 

r=l r=l s^r 

X0 (-q„ • gr rs ) (q rs ■ g rs f q rs q rs / d-y 5 (r - q r + 7<? rs ) . 

Jo 



(Dll) 



The first term on the right hand side is the called kinetic part of the momentum flux, while 
the second term is the collisional transfer part. 
Finally, the energy density is given by 



N 



SiV-r) = Y J ^(r-qr) 



(D12) 



Proceeding in a similar way as for the momentum density, 

N 2 N N 

L£ (r; r) = ^v r 8(r - q r ) + y ]T ]T T 1 (r, s) [t# (r - q r ) + 7; s 2 5 (r - q.j\ 
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Emv r m(l + a)cr v ^^ _ 
-^-v r 5(r - q r ) + 2^ 2^ 5 ^ rs ~ a >® (~ qrs ' 9rs > 



r=l 



r=l s^r 



X (Qrs ■ g rS ) 2 Qrs ■ GrsQrs I S (r - q r + ^q rs ) 

Jo 



777(1 — a 2 ) 



N N 



XI ~ o-)0 {-q rs ■ g rs ) \q rs ■ 9rs\ 3 S (r - q r ) . (D13) 



This consists of two parts, one which can be expressed as a gradient and a second one that is 
inherently local, characterizing the loss in energy due to the inelastic character of collisions. 
The energy balance equation becomes 



J^(r ; r, f) + V • 8 (r ; r, t) = -w (r ; r, t) , 



(D14) 



with the energy flux density 



A? 



r=l 



N N 



(r;r) = X,~9 -M(r-g r )+ ^^(5(g rs -rr) 



r=l s^r 



x6 (-q„ • gij) (q rs ■ g rs ) 2 (q rs ■ G rs ) q rs [ d'y 5 (r - q r + ^q rs ) , (D15) 

Jo 
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and the energy source 

r=l s^r 

This completes the identification of the flux and source terms in the microscopic balance 
equations for the forward dynamics. 

A similar analysis can be done for the balance equations associated with the backward 
dynamics of the phase functions, described by the equations 

dAa d(-t) ^ " L Aa (F; r ' ~ t] = °' (D17) 
where t > 0. The form of the conservation laws now is 

dAa d{-t) ^ + V ' j « (F; r ' ~ t] = W « (F; r ' ~ t] ■ (D18) 
If the generators L and L_ were the same, then this would just be the time reversal of 
Eq. (D2). However, the contributions for positive times from the binary collision operators 
T are replaced by those from T_ when considering negative times. Repeating the 
above analysis for this case, gives the time reversed fluxes as 

n m(l + a)a N N 
h~(r;r) = m ^ v r v r 5 (r — q r ) H — ^^5(g rs -cx) 

r=l r=l s^r 

xQ(q rs - g rs )(q rs - g rs ) 2 q rs q rs [ d'y 5{r - q r + 7<? rs ), (D19) 

Jo 

N , s N N 



r=l r=l s^r 

Qrs ) \Qrs 

■ G rs )q rs / d^5(r - q r + 

Jo 



(D20) 



while the expression for the reverse energy source is 
m (l _ a z\ N N 

w~(T;r) = — ^2^2S(q rs - a)Q(q rs ■ g rs )(q rs ■ g rs ) 3 S (r - q r ) . (D21) 

r=l Sy^r 

The superscript — is used here to indicate that the quantities are associated with the time 
reversed equations. The source term w~ in the balance equation for the energy is now 
positive, accounting for collisional increase in energy on moving backwards along a trajectory 
originally generated forward in time. In the elastic limit, this latter effect vanishes. However, 
the fluxes still differ from those for t > 0, with the collisional transfer contributions being 
defined on different pre-collision hemispheres. 
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1. Dimensionless Balance Equations 



The purpose of this subsection is two fold. First, the dimensionless forms of the vari- 
ous fluxes that determine the hydrodynamic parameters are identified. Second, a special 
property of the source term in the equation for the energy noted in ref. [3], namely that 
it is orthogonal to the set of conjugate densities m the homogeneous limit, is shortly 

reviewed here. 

Dimensionless forms for the microscopic densities Af, Q, and £ are defined by (recall the 
choice / = n^ 1 ^ made for the length scale) 

^ (r> ,,^. £>._,,), (D22) 

£*(r*; r*) = = f>f5(r* - q* r ), (D23) 



n h T h j 
0(r;r) N 



0*(r'; r*) = Vl ' ' = £ " # )■ 

mn h v {T h ) 

With this choice, the balance equations in the dimensionless form become 



(D24) 



dN* d 

-S-+8?- <r = ' (D25) 
|-f) e ' + 4' h "°' <D26) 

I - c ») £ ' + 7^ - s ' = - w '- (D27) 

The functional forms of the fluxes and source terms are related to those given above by 

h "( r *-*)-^Pu^ |h(r; '- )l -'--' (D28) 

u>-(r«; r«) = ^ T ;jl = 2 [«,(r; r)] r , r .^. . (D30) 

It is understood that, in addition to the indicated substitutions, the changes a — > er* = a /I 
and m — > 1, are made. The direct densities functions a* a (T*;r*) appearing in the response 
functions are defined by Eq. (17). In the case of a* 2 , the balance equations give 

-~r E \-r-8 = — t- ( D31 ) 



ds d dr* \ d J d 
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For reasons that become apparent below, add a term J2 a W a *a 011 both sides of the 
above equation to write 



ds 



where the new source £* is 



dr* 



d 



e*(F*;r*,s) 



2w* 
~cT 



-Co 



3 * ( <91nCo . .. 



Evidently, the Fourier representation of Eq. (D32) is 



|i + k£* (o) a a - ik *.( 2 -?--g*)= t (fc*, s) 



d 



(D32) 



(D33) 



(D34) 



The homogeneous limit of the source term, l*(Q,s), has an important orthogonality 
property, which is the reason for adding the contribution J2 a W a *a above. To see 
this, the Fourier transform for k* = of the second term on the right hand side of Eq. 
(D33) is written in the equivalent form 



= e 



sC* 



^«(o) + Co(|^ + i)^(o) 



(D35) 



The last equality is obtained as follows. From the definition of ip* in Eq. (19) it is 
2 



V*d 



I 



dr* w* (o) c (o) 



^ r; c-n \i-Nd N a l 
d l ' V ° (Th)] n h T hM T k )V 



dTw(T;0) 



x 



dp h (T;n h ,T h ) dp{T;n h ,T h ) 

-Oal H 7^ 0a2 



dn h 



dT h 
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n h T h v (T h )Vd \ 6ainh dn h 0< " lh dT h 



x 



J dTw(T;0)p h (T;n h ,T h ). 



Moreover, from Eq. (D14), 

/d 
drw(T;0)p h (T;n h ,T h ) = -n h V( (n h ,T h )T h , 

so that Eq. (D36) can be rewritten as 

d\n( (n h ,T h ) 



V*d 



I 



dr w* (o) c (o) = a Sai 



1 + 



d In rih 



+ Sa2~ 



(D36) 



(D37) 



(D38) 
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This proves Eq. (D35). Use of it into the Fourier transform of Eq. (D33) particularized for 
k* = yields 

? (r*; 0, s) = e sC * (1 - P f ) ^w* (0) (D39) 
with the operator pt given by 

p f x (r*) = v*- 1 ( r *; °) / dT * C (r*; o) x (r*) . (D40) 

This is a projection operator since the quantities {a*(r*;0)} and |^(r*;0)| form a 
biorthogonal set, in the sense that 

V*- 1 J dT*Z* a (T*;0)^(T*;0) = 6 a(3 , (D41) 

as shown in ref. [3]. The utility of this result is that t (0, s) is the source phase function 
for the transport coefficient ( u , as shown in Eq. (F2). Then, the presence of this orthogonal 
projection is essential for the existence of the large s limit in the Green-Kubo representation 
for ( u , as discussed in reference [3]. A similar orthogonal projection occurs for all of the 
direct fluxes F s '- f occurring in the representation (32) for transport coefficients. 

In summary, the conservation laws associated with the chosen phase functions {a a } in 
dimensionless variables take the form 



cfn* (V*- k* si ~ ~ 

r ' +EC ( r " °) ( r " s ) ~ ik * ■ ^ (r*; fc*, s) = s a2 t (r* ; k\ s ) , (D42) 
as p 

where K*^ d (T*; 0) is the k* = limit of the matrix given in Eqs. (11)-(13) and 

f*(T*;r*) = g*(T*;r*), (D43) 

/*(r*; r*) = ^s*(T*; r*) - g*(T*; r*), (D44) 

f 3 ^.(r*;r*) = h* J (r*;r*). (D45) 

The last equation above gives the tensor flux appearing in the equation for the dimensionless 
momentum Q*{Y*- r*) = a* 3 (F*; r*). 

APPENDIX E: CONJUGATE FUNCTIONS AND CONSERVATION LAWS 

Consider the conjugate functions ^* (k*) that occur in the response function (15), ob- 
tained from the linear response analysis, and are defined by Eq. (19). The dimensionless 
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dynamical equation obeyed by these functions is 

(d s +r)r a (k*, S ) = o. (ei) 

In the homogeneous limit k* = 0, these functions verify [3] 

rc(o) = ^/c;r(o)^(o), (E2) 

where K.* hyd (0) is the hydrodynamic transport matrix given by Eqs. (11) and (13) in the 
main text. This property allows the construction of a new set of "conservation laws": 

|-c (fc*, s) + J2 C (°) ^ ( fc *> s ) - • ^ ( fc *> s ) = °' ( E3 ) 

OS p 

with the "conjugate fluxes" 7* (fe, s) defined by 

l k* ■ 7; (fc', s) = - (r - £ /C?/ (0) j ^ (fc*, s) . (E4) 

The property (E2) assures that ifc* • 7* vanishes at fc* = 0, justifying the introduction of 
7* as a flux. It is the conjugate quantity occurring in the Green-Kubo expressions of the 
transport coefficients, Eq. (35). More precisely, T* in that expression is defined by 

T; = k* ■ 7; (0) . (E5) 



APPENDIX F: TRANSPORT COEFFICIENTS: SOME DETAILS 

The linear response analysis of ref. [3], identified the transport coefficients in the inter- 
mediate Helfand representation by Eqs. (118)-(120). The transformation of those results to 
the dimensionless, hard sphere form is outlined in Sec. VIII of that reference. Some further 
details for the expressions of the dimensionless transport coefficients are described briefly 
here. 

The dimensionless response functions have a dynamics generated by C , as in Eq. (B18). 
However, the transport coefficients are given in terms of correlation functions whose dy- 
namics is generated by C* — \jC* hyd (0)] , where the superscript T indicates transposed. As 
discussed in ref. [3] and also in the main text here, this generator compensates for both 
the cooling of the reference state (through the change of L* by C*) and the dynamics of 
homogeneous perturbations of that state (through the subtraction of K,* hyd (0)). 
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Consider first the Euler transport coefficient ( u given in the intermediate Helfand repre- 
sentation by Eq. (118) of ref. [3], whose dimensionless form is 

C^ = -limfc*.^ 1) ( a ), (Fl) 

where S 23 (s) is obtained from 

sl 3 (k*,s) = llf^l ^ 1 J dr*t(r*,k*)e- sT ^(r*;-k*) 

= V*- 1 J dT* I* (r*; k*) e~ S ( £ (r*; -Jfe*) , (F2) 

as explained below. In this expression, £* (T*; k*) is the source term in the balance equation 
for the energy given in Eq. (D33) and ip*(T*; k*) is the conjugate density associated with 
the longitudinal component of the flow velocity. Equation (F2) is easily derived by using 
Eq. (168) in ref. [3]. It can also be expressed as 

S* 23 {k*, s) = V*~ l J dT* t{T*, k*)$*{T*; -k*, s), (F3) 

with 

^|(r*; k*, s) = ^^ Q ( S )C(r ; k*), (F4) 



the matrix evolution operator li*(s) being given by 



W = Cff-\-^ = e^-"W e — *. (F5) 

J pa 

Here, K,* hyd (k*) is as always the matrix given in Eqs. (11)-(13). It will be seen in the following 
that, consistently with the description given at the beginning of this Section, going to the 
dimensionless expressions of the transport coefficients implies replacing the matrix evolution 
U(t,T) used in ref. [3] by U*(s) defined above. 

The quantity S* 23 \s) in Eq. (Fl) is generated from S* 23 (k*,s) by means of the Taylor 
series expansion around k* = 0, that for any function X(k*) is defined as 

X(k*) = X(0) + ik* ■ X™ - k*k* : X< 2 > + . . . . (F6) 

Then, it is 



,sK.* hyd (0) 



$23 



(1) 



(s) = -V*- 1 J dT*£*{T*;0)e + < °) ^(T*), (F7) 
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with 

N r, 

^(l) (r) = _^ q; __ p * (r) . (F8) 

r=l S M 

Upon writing Eq. (F7), a contribution involving £*^(T*) has been omitted since it vanishes 
due to the symmetry of ip*(T*;0). Taking now into account Eq. (D39) and the definitions 
given in (D30) and (38), it is easily verified that 

r(r*;o) = (i -p^)w*(t*). (F9) 

Use of this into Eq. (F7) and later substitution of the result into Eq. (Fl) leads to Eqs. 
(47)- (50), after using the symmetry of the integrand. 

The Navier-Stokes transport coefficients at order k 2 associated with the heat and mo- 
mentum flux, are given by Eq. (119) of ref. [3] and, more explicitly, they are identified in 
Appendix F of that reference. In dimensionless form, they are the elements of the matrix 

A* = lim k*k* : [D* (1) (s) - D*(0, 0) C (1> (s)l , (F10) 

where the correlation functions C*(k*, s) and D (k*, s) are defined by 

C* ap (k* , a) = V*- 1 J dT*a* a (T*;k*)^(T*;-k*,s). (Fll) 

D* ap (k*, s) = V*- 1 J dT* / * (r*; k*) (r ; -k*, s) , (F12) 

with the direct densities a* given in Eq. (17) and the associated fluxes f* given in Eqs. 
(D43)-(D45). 

The explicit form of the first order in k* expanded correlation matrix C (s) above is 

c*S\s) = v*- 1 J dr*a*U(r*)r p (r*;o, s ) 

-V*' 1 J dT*a* a (r;0)^ {1 \r;s) 
= -V*- 1 J dT* a* a (r*; 0) (r*; s) . (F13) 

The contribution involving a^^T*) is proportional to the center of mass coordinates and 
vanishes from symmetry. A similar analysis of D**^ (s) gives 

< 1} (*) = -v*' 1 J dr* r a (r ; o) (r ; 8 ) . (fm) 
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The two terms on the right side of Eq. (F10) now combine and the transport matrix A* 
takes the simple form 

A* af3 = - limfc*fc* : V*- 1 J dT* /*(P; 0)(1 - P)^J (1) (P; s) , (F15) 

where P is the projection operator 

PX(T*) = C( r *; O)^ 1 / dT*a* a (T*; 0)X{T*). (F16) 

a J 

Since the operator 1 — P is orthogonal to the invariants of U*(s), as shown in ref. [3], Eq. 
(F15) is equivalent to 

A^ = -limV*- 1 J dT*F* f (T*)k* ■ ^J (1) (P; s), (F17) 

where F*f is the projected flux 

F *f(F*) = (1 - Pt) %* ■ f* (r*; 0) . (F18) 

The projection operator P 1 " is the adjoint of P, and it is the same as that given by Eq. 
(D40). 

1. Shear Viscosity 

The dimensionless shear viscosity is determined from A* 4A through 

v * = A* 4 = -UmV^Pk* : J dT*F* f (T*) U* (P) ^ 

= -limr- 1 fefe: y"rfr*^(r*;0))e" S ( £ + ^)^ 4 * (1) (r*). (F19) 

The projection operator in the definition of the projected fluxes, Eq. (F18), does not con- 
tribute in this case due to the symmetry of the integrand. The conjugate moment, ipl^ (P), 
is easily obtained from the definition in Eq. (F6) and Eq. (25), 

k* ■ (r*) = - £ (k* ■ q ;) -£- P \ (r*) , (F20) 

where v* ±1 is the first of the transversal components of v r . Moreover, using Eq. (D45), 



fe-/;(r*;0) = Hi ±1 (r), (F21) 
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so Eq. (F19) becomes 

f-^fvKjr^W), (F22) 

with 

r=l r 'f 

The x and y components in the above expression are due to the arbitrary choice of the x 
axis along k and the y axis along the first tranverse direction. 

An alternative expression for the shear viscosity can be derived as follows. Define 

M^n = -\ £ + - \k 3 € • ^) P*(n, (F24) 

where i and j are spatial coordinates. The symmetry of Eq. (F22) implies that it is equivalent 
to 

= - d 2 + 1 d+2 ^ v *~ 1 J2J2 f d^une~ s ^ + ^M;^n. (f 25 ) 

i=i j=i 

This is the expression presented in Eqs. (56) and (57) in the main text. 



2. Bulk Viscosity 

The bulk viscosity n* occurs in the 33 matrix element of the hydrodynamic transport 
matrix, namely it is 

K * + 2(yd ~ l)r] * =a* 33 = - limy*- 1 J rfr*F 3 *V s ( £ * + ^)^*(r). (F26) 

The projected flux is found to be 

F/(r*) = k* J; (r*5 0) - ai (r ; 0) - f a 2 (r ; 0) (F27) 

z c m z 

and the conjugate moment is 

A<5(r') = fc* • v 3 * (1) (r*) = - E (fc* • tf) fc* • ^ph (n ■ (F28) 

r=l r 

Taking into account that k* ■ / 3 * (r*; 0) = H*, 3 (r*) and the symmetry of the tensor H*, the 
expressions used in Sec. Ill, Eqs. (68) and (69), follow directly. 
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3. Thermal conductivity 



The dimensionless thermal conductivity A* is given by 



d 



d 



A* = ^A* 22 = -^MmV*- 1 I dT*F* f 



u{s)k* • v* (1) (r*) 



J 2 



( r _io 



--limV*- 1 / ■/! / :; ; r 



fc • v 2 * (1) (r*) 



The projected flux is 



F? = 



k* ■ 



with 



f*(F*;0)--^ P* 



s*(T*;0)-[p* h + ^\P* 



N 



r=l 



and the conjugate moment is 

k* ■ (r) = -\Y. (** • <?*) ^ • [«(r) 

r=l r 

The above expressions are equivalent to Eqs. (75) and (76). 



N 



(F29) 



(F30) 



(F31) 



(F32) 



4. The /j, Coefficient 



The expression of the coefficient ji* from linear response reads 
- d 



-A* 

, iv 21 



= —Hmv*- 1 J dr*F; f [w*(s)fc* • v* (1) (r*) 



= _^i imV — 1 / rfr*F;V sr 



_ ;j! r , o 9 In Co-* _ ^*(i)(p*) 



fc*-<0* (1) (r*)-2 



—d— — — lim V 



<91nn^ 

T> k* • v5 (1) (r* 



(F33) 



Slnn/j 

Since the projected flux is the same as in Eq. (F29), the second term on the right hand side 
of Eq. (F33) is proportional to A*, and the expression can be written 
3 In Co 



JT = fj, — 2 
d 



d In rih 



A* 



7 



--limy*- 1 / dr*F; f e - sC k* 



(T*) - 2^-k* ■ 



(F34) 
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Therefore, the involved conjugate moments are 



k* ■ ii>l {1) (r* ) = [lv (T)] iva n I drk* r 



Nd 



k* • v 2 * (1) (r*) 



[lv (T)] Nd T J 



dr k* ■ r* 



'Spih(T\ {y a }) 
5n(r) 

Spi h (T\ {y a }) 
5T(r) 



{y a }={n,T,0} 



{y a }={n,T,0} 



(F35) 



(F36) 



Substitution of Eqs. (F35) and (F36) into Eq. (F34), and use of the symmetry of the inte- 
grand leads to Eqs. (84)-(86). 



APPENDIX G: ELASTIC HARD SPHERES 



The aim of this Appendix is to give the Green-Kubo expressions for the transport co- 
efficients in the case of elastic hard spheres or disks. The linear response method used in 
reference [3], was specifically applied to initial perturbations of a local HCS. The results ap- 
ply as well in the elastic case for perturbations of a local equilibrium distribution, with only 
the replacement of the local HCS ensemble by a local equilibrium Gibbs ensemble. Here, for 
simplicity and to compare with standard results in the literature, a local grand canonical 
ensemble will be used. Strictly speaking, this is not the a — > 1 limit of the results derived 
in this paper for a granular fluid, since the HCS ensemble in the elastic limit goes over to a 
microcanonical ensemble with constant total energy and momentum, i.e., the so-called MD 
ensemble. The results obtained by using different ensembles are known to be equivalent up 
to fluctuations that do not contribute in the thermodynamic limit, but the ensemble used 
affects the forms of the projected fluxes appearing in the Green-Kubo expressions. 

The form of the representations for all transport coefficients is, of course, the same as in 
the inelastic case, but with some simplifications. The most obvious ones are the absence of 
cooling, implying that C* reduces to L , and the absence of any dynamics for homogeneous 
perturbations, that reflects itself in that K,* hyd (0) = 0. Of course, both simplifications are 
closely related. Moreover, the dimensionless time s is now simply proportional to t. The 
direct fluxes F* s in Eq. (35) vanish since all the the source terms are zero and, consequently, 
also all the transport coefficients related to the cooling rate. On the other hand, the direct 
fluxes F*f are unchanged. The generic transport coefficient in Eq. (90) then becomes 

u e i = limV*- 1 [ dY* F f (P) e- sT M* (P) . (Gl) 
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It remains only to determine the moments M* (T*), and afterwards the Green-Kubo rep- 
resentation can be constructed directly as done in the main text for the general inelastic 
case. 

The distribution of the local grand canonical ensemble is [13], 

s(r-r) 



p lGC (T)=Z(T)expt-Q L + / dr 



g(r)Af(T;r) 



T(r) 



+ Y(r)-g(T;r) 



, (G2) 



where the parameters of the ensemble are defined as 



g(r) = 



T{r) 
Y{r) = 



?(r)- 



mU 2 (r) 



2T(r) 

n(r)mU(r) 
T(r) ' 



(G3) 



(G4) 



q being the chemical potential, U the flow velocity, T the temperature in units where the 
Boltzmann constant is unity, and Ql the normalization constant. Moreover, S(r) is the 
overlap function that is zero for any configuration of overlapping pairs and unity otherwise. 

Consider first the conjugate fluxes defined in Eq. (19) with the substitution of pih by 
Pigc, i- e., 



r a (T*;k*) = M a / dre 



Akr 



Spice [£j {y/s}] 
$y a (r) 



(G5) 



{yp}={n,T,0} 



The functional differentiations can be performed explicitly in this case, with the results 



^(r*;fc*)= l + n h h e (k) Af*(T*;k*)p* GC (T*), 



(G6) 



d 



r 2 (T*;k*)= £*(T*;k*)--jV*(r*;k*) p GC (r*), (G7) 

^ (r*; k*) = 2g(T*; k*)p GC (T*). (G8) 
Here h e (k) is the Fourier transform of g e (f) — 1, where (7 e (r) is the equilibrium radial 



distribution function. The notation if)^ = <^tppif}* ± j has been used. It is seen that the 
■ip a (r*; fe)'s are linear combinations of the direct fluxes a a (T*; k) multiplied by the grand 
canonical ensemble dimensionless distribution p* GC (T*). 

The moments M* (T*) appearing in the Helfand expressions, are the coefficients of ik* in 
the expansions of the ip a (fe*)'s, as indicated in Eq. (33), 



M*(T*) =k* ■ V* (1) (r*). 



(G9) 
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In detail, they are found to be 

N 



p.^f\r)=[l + n h h e (0)y 1 J2k*.q;.p* GC (r), (G10) 

r=l 

k* ■ v 2 * (1) (n = E** • 9; < 2 - 2 ^ c (r*), (en) 

r=l ^ ' 

fc* • ^; {1) (r) = 2^ (** • <?;) (** • <) Peer), (G12) 



r=l 

AT 



fc* • i>$\n = 2j2(fr- <?;) <^c(n- (gw) 

r=l 

Once these moments have been identified, the various transport coefficients for a system of 
elastic hard spheres or disks can be directly written. 

1. Shear and Bulk Viscosities 

The Helfand forms in Sec. Ill for the shear and bulk viscosities of the granular model, in 
the elastic limit become 

rj* = lim Q V H (s) , k* = lim Q K H (s) , (G14) 
respectively, where Q n H (s) and Vt^ (s) are now equilibrium time correlation functions, 

(«) = - d Z*d_2 ^^j dr (KaPoc) . (G15) 



i=i i=i 



fi^ ( s ) = -2(\/*rf 2 )^ 1 y rfr* tr H^e"' 1 -* (M*p* GC ) , (G16) 
with given by Eq. (43) and H?- by Eq. (71). Moreover, 

N / 2 \ 

K« = E kx* + tfx* - 2 5 «* • v * ) ' (G17) 

r=l ^ ' 



N 



K = (Gi8) 



These are the usual well-known results for systems of elastic hard spheres or disks, in di- 
mensionless form. The corresponding Green-Kubo expressions are 

rj* = Vl v + lim f ds'n r G (s') , (G19) 
Jo 
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k* — + lim / ds'tt K G (s') . (G20) 
Jo 

The instantaneous parts, Q$ = Q V H (0) and Qq = f2^(0), now can be evaluated exactly with 
the result 

"S = — "3= W * (g) - (G21) 

The calculation of the conjugate fluxes in the expressions of Q G (s) and Q G (s) requires some 
care. The calculation of T* ■„■ has been discussed in detail in Sec. Ill, and the result is 
given by Eq. (63). The subsequent expression for the elastic shear viscosity is presented in 
Eq. (65). In a completely similar way, it is obtained that the integrand in the Green-Kubo 
expression for the bulk viscosity is 

tt G (s) = {V*d 2 )~ l J dT* tr H* / (r*)e" sr [tr H*-(P>^ c (r*)] . (G22) 
2. Thermal Conductivity 

The expression of the thermal conductivity in the Helhand representation reads 

A* = lim^(s), (G23) 

where for a system of elastic hard spheres or disks, is the time equilibrium correlation 
function 

n x H (s) = -^d)' 1 J dT* S* f ■ e- sT (M* x p* GC ) . (G24) 
The moment M{ is identified from Eq. (Gil) as 

X ' d 



M A * = J>* < 2 -^)- (G25) 

r=l ^ 

The corresponding Green-Kubo representation is 



A* = SI* + lim f ds' Vl G (s') . (G26) 
Jo 

Again, the instantaneous contribution Oq = f2#(0) can be calculated exactly, 

n (d-l)/2 a *d+l 

The analysis of the Green-Kubo flux Ta is similar to that for the viscosities, leading to the 
result 

tt G (s) = {V^d)- 1 J dT* S* f ■ e- sV (S*-p* GC ) , (G27) 
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where S*~(T*) is the volume integrated energy flux for the time reversed microscopic laws 
identified in Appendix D as 

N 

S*-(T*) = J2< 2V * 
r=l 

N N 

+(T * E E 5 ^s ' (Q*rs ■ 9*rs) (Q*rs ' 9*rsf (Frs ' GJJfc. (G28) 

r=l s^r 



3. /x Coefficient 

It has been mentioned several time that the transport coefficient ji vanishes in the elastic 
case, and it is instructive to see this at the formally exact level of linear response theory. 
The elastic limit of the Helfand form for /x* in Eq. (F33) is 



where 



~ lim V*- 1 J dT* S* f ■ e~ sL * (m; P * gc ) , (G29) 



-l N 

M; = ^; (1) = l + n h h(0)\ ^g r *. (G30) 



r=l 



This is proportional to the center of mass position of the system. Its time evolution is, 
therefore, proportional to the total momentum P*. But, by definition, S* f is projected 
orthogonal to P*. As, moreover, it is easily seen that the average on the right hand side of 
Eq. (G29) vanishes for s = 0, it is concluded that it vanishes for all times. 



APPENDIX H: STATIONARY REPRESENTATION FOR MD SIMULATIONS 

The transformation to dimensionless velocities defined in the text, is based on scaling 
relative to v (t) = [2T h (t)/m] 1/2 . This was done so as to be able to pose theoretical questions 
of interest in an elegant self consistent form. However, this is inconvenient in practice since 
the cooling rate is given implicitly in terms of the stationary HCS. Instead, the same analysis 
can be performed by scaling with a known function £(t) instead of v (t) [5], to get a Liouville 
equation in the form 

J^p**(r**, t) - dln d f^ E ^ ' [«*V*(r**, *)] + ^i**(r*>**(r**, *) = o, (hi) 

r=l 
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where T** = {q**,v**} = v r /£(i)} • Here, as earlier, / is a constant characteristic 

length in the system, 

P »(r~,t) = [e x (t)] Nd p(r,t), (m) 

and 

r*(r*) = ^I(r) = [I(r)] Ivr ... (h 3 ) 

Next define a new time variable r by 

dr = (H4) 

and choose £(t) to make the coefficients of the Liouville equation (HI) independent of r, 
namely verifying 

where £o is an arbitrary dimensionless constant that can be picked for convenience. Then, 
the choice = 2l/£ t is made and Eq. (H4) becomes 

*=-. (H6) 

The scaled Liouville equation (HI) takes the form 

3 P N a 

Q-/* + | E ■ «P**) + L *V* = 0. (H7) 

r=l r 

This is formally the same as Eq. (B2), except that here the cooling rate has been replaced 
by the arbitrary constant £ - 

There is a stationary solution p** t to Eq. (H7) determined by 



C (r";&)pS(r";&)=0, (H8) 

r-(r-? ) = f£ A. „;+!•-. (H9) 

r=l r 

Clearly p**(r**;£ ) is the same function as p* h (T*) with only the unknown value Q replaced 
by the constant £ • Interestingly, it is possible to determine Co from the chosen value of £ 
and the measured value of the temperature T* t * of the steady state, defined by 

T* t * = i J dV** v{* 2 p** (r**; f ). (H10) 
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The three mentioned quantities are related by 



Co 



(Hll) 



This relationship may be derived as follows. Define for homogeneous systems in general the 
temperature 

T{t) 



T**(r) = ^J rfr**<V*( r **> r ^o) 



(H12) 

m X 2 {t) 

i.e., the temperature T** is the actual temperature but expressed in the arbitrary scaling 
variables. Now consider the dynamical equation associated with this scaled temperature, 
that follows directly from Eq. (H10), 



d 



-Co )T**(t) 



Ko(t)T**(r) 



(H13) 



dr "7 * > m ' 

Using now the scaling property of the distribution function of the HCS and, therefore, of 
p**(r**;£ ) this becomes 



(H14) 



whose solution is 



T**(r) 



f 2 

2Co 2 



r?o/2 



(H15) 



CoV 2T **(°)- 

which in the long time limit goes to Eq. (Hll) above. Therefore, in practice, one imagines 
measuring T s * t * rather than solving for Q self-consistently in the HCS state. Also, the 
different generators defined earlier and the stationary representation of two-time correlation 
functions over the HCS ensemble can be translated into this language of arbitrary scaling 
[5, 6]. 
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